Advertisement
Guest User

Untitled

a guest
Jul 27th, 2017
60
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 5.87 KB | None | 0 0
  1. // orbit - Program to compute the orbit of a comet.
  2.  
  3.  
  4. /***************************************************************
  5.                 -- Code to use gnuplot --
  6. Include these lines of code at the top of your c or c++ program
  7. to use gnuplot for graphical output.  Call open_gnupipe(), then
  8. use gnuplot("command string") to generate plots.  Finally call
  9. close_gnupipe() before exiting the program.  This routine
  10. assumes that your path includes the gnuplot program.
  11.                     (Dr. Mark C. Fair)
  12. ***************************************************************/
  13. #include <stdlib.h>
  14. #include <stdio.h>
  15. #include <string.h>
  16. FILE *gnu_pipe;
  17. void open_gnupipe(){
  18.     if (( gnu_pipe = popen("gnuplot -persist","w")) == NULL){
  19.         perror("popen failed");
  20.         exit(1); }
  21. }
  22. void gnuplot(char *cmdstr){
  23.     char cmdstr1[strlen(cmdstr)+1];
  24.     sprintf(cmdstr1,"%s\n",cmdstr);
  25.     fputs(cmdstr1,gnu_pipe);
  26.     fflush(gnu_pipe);
  27. }
  28. void close_gnupipe(){
  29.     if (pclose(gnu_pipe) == -1)
  30.         printf("Problem closing communication to Gnuplot\n");
  31. }
  32. /***************************************************************
  33.                 -- End of code to use gnuplot --
  34. ***************************************************************/
  35.  
  36.  
  37.  #include <iostream>
  38.  #include <fstream>
  39.  #include <cmath>
  40.  #include <cassert>
  41.  
  42.  
  43.  using std::cout;
  44.  using std::cin;
  45.  using std::endl;
  46.  using std::ofstream;
  47.  
  48.  
  49. void gravrk(double*,double,double*,double*);
  50. void rk4( double*,int,double,double,
  51.     void (*derivsRK)(double*,double,double*,double*),
  52.     double*);
  53.    
  54.  
  55.  int main()
  56.  {
  57.  
  58.     double param[5];
  59.     param[0] = 0.5;
  60.     param[1] = 1;
  61.     param[2] = .05;
  62.     param[3] = .1;
  63.     param[4] = .01;
  64.    
  65.     double state[4];
  66.     double nState = 4;
  67.     double tau = .1;
  68.     double time = 0;
  69.     int nStep = 400;
  70.     double *tplot,*rplot,*thplot;
  71.     tplot = new double[nStep];
  72.     rplot = new double[nStep];
  73.     thplot = new double[nStep];
  74.    
  75.    
  76.     double r[2],v[2];
  77.     v[0] = 0;
  78.     v[1] = 0;
  79.    
  80.    
  81.     cout << "enter in the initial height(cm): ";
  82.     cin >> r[0];
  83.     cout << "enter in the initial angle(rads): ";
  84.     cin >> r[1];
  85.    
  86.     state[0] = r[0];
  87.     state[1] = r[1];
  88.     state[2] = v[0];
  89.     state[3] = v[1];
  90.    
  91.     for(int iStep = 0; iStep < nStep; ++iStep){
  92.    
  93.             rplot[iStep] = r[0];               // Record position for plotting
  94.             thplot[iStep] = r[1];
  95.         tplot[iStep] = time;
  96.    
  97.         rk4( state, nState, time, tau, gravrk, param );
  98.         r[0] = state[0]; r[1] = state[1];   // 4th order Runge-Kutta
  99.         v[0] = state[2]; v[1] = state[3];      
  100.         time += tau;
  101.    
  102.     }
  103.    
  104.     ofstream fout("12.txt");
  105.    
  106.     for(int i = 0; i < nStep; ++i)
  107.     {
  108.         fout << rplot[i] << " " << thplot[i] << " " << tplot[i] << endl;
  109.    
  110.     }
  111.    
  112.    
  113.     open_gnupipe();
  114.     gnuplot((char *)"set term x11 3 font 'Times,14'; \
  115.             set title 'wilberforce pendulum' font 'Times,18'; \
  116.             set grid; \
  117.             set xlabel 'Time(s)' font 'Times,18'; \
  118.             set ylabel 'z(cm) theta(radians)' font 'Times,18'; \
  119.             plot '12.txt' using 3:1 with lines title 'Z(t)'; \
  120.             replot '12.txt' using 3:2 with lines title 'theta(t)';\
  121.             set term postscript eps color 'Times' 18; \
  122.             set output '12.eps'; \
  123.             replot");  
  124.     close_gnupipe();
  125.    
  126.    
  127.     delete [] tplot,rplot,thplot;
  128.    
  129.     return 0;
  130.  }
  131.  
  132.      
  133. void gravrk(double x[], double t, double param[], double deriv[]) {
  134. //  Returns right-hand side of Kepler ODE; used by Runge-Kutta routines
  135. //  Inputs
  136. //    x      State vector [r(1) r(2) v(1) v(2)]
  137. //    t      Time (not used)
  138. //    param     mass, I, k, delta, e
  139. //  Output
  140. //    deriv  Derivatives [dr(1)/dt dr(2)/dt dv(1)/dt dv(2)/dt]
  141.  
  142.   //* Compute acceleration
  143.   double mass = param[0];
  144.   double I = param[1];
  145.   double k = param[2];
  146.   double delta = param[3];
  147.   double e = param[4];
  148.  
  149.   double r1 = x[0], r2 = x[1];     // Unravel the vector s into
  150.   double v1 = x[2], v2 = x[3];     // position and velocity
  151.  
  152.   double accelz = -(k/mass)*r1 - .5*(1.0/mass)*e*r2;  // Gravitational acceleration
  153.   double acceltheta = -(delta/I)*r2 - .5*(1.0/I)*e*r1;  
  154.  
  155.   cout << "v1: " << v1 << endl; //<< "v2: " << v2 << endl;
  156.   //cout << "accelz: " << accelz << endl << "acceltheta: " << acceltheta << endl;
  157.  
  158.   //* Return derivatives [dr[1]/dt dr[2]/dt dv[1]/dt dv[2]/dt]
  159.   deriv[0] = v1;       deriv[2] = v2;
  160.   deriv[1] = accelz;   deriv[3] = acceltheta;
  161. }
  162.  
  163. void rk4(double x[], int nX, double t, double tau,
  164.   void (*derivsRK)(double x[], double t, double param[], double deriv[]),
  165.   double param[]) {
  166. // Runge-Kutta integrator (4th order)
  167. // Inputs
  168. //   x          Current value of dependent variable
  169. //   nX         Number of elements in dependent variable x
  170. //   t          Independent variable (usually time)
  171. //   tau        Step size (usually time step)
  172. //   derivsRK   Right hand side of the ODE; derivsRK is the
  173. //              name of the function which returns dx/dt
  174. //              Calling format derivsRK(x,t,param,dxdt).
  175. //   param      Extra parameters passed to derivsRK
  176. // Output
  177. //   x          New value of x after a step of size tau
  178.  
  179.   double *F1, *F2, *F3, *F4, *xtemp;
  180.   F1 = new double [nX];  F2 = new double [nX];
  181.   F3 = new double [nX];  F4 = new double [nX];
  182.   xtemp = new double [nX];
  183.  
  184.   //* Evaluate F1 = f(x,t).
  185.   (*derivsRK)( x, t, param, F1 );  
  186.  
  187.   //* Evaluate F2 = f( x+tau*F1/2, t+tau/2 ).
  188.   double half_tau = 0.5*tau;
  189.   double t_half = t + half_tau;
  190.  
  191.   for(int i=0; i< nX; ++i)
  192.     xtemp[i] = x[i] + half_tau*F1[i];
  193.   (*derivsRK)( xtemp, t_half, param, F2 );  
  194.  
  195.   //* Evaluate F3 = f( x+tau*F2/2, t+tau/2 ).
  196.   for(int i=0; i< nX; ++i)
  197.     xtemp[i] = x[i] + half_tau*F2[i];
  198.   (*derivsRK)( xtemp, t_half, param, F3 );
  199.  
  200.   //* Evaluate F4 = f( x+tau*F3, t+tau ).
  201.   double t_full = t + tau;
  202.   for(int i=0; i< nX; ++i)
  203.     xtemp[i] = x[i] + tau*F3[i];
  204.   (*derivsRK)( xtemp, t_full, param, F4 );
  205.  
  206.   //* Return x(t+tau) computed from fourth-order R-K.
  207.   for(int i=0; i< nX; ++i)
  208.     x[i] += tau/6.*(F1[i] + F4[i] + 2.*(F2[i]+F3[i]));
  209.    
  210.   delete [] F1, F2, F3, F4, xtemp;
  211. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement