Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- // orbit - Program to compute the orbit of a comet.
- /***************************************************************
- -- Code to use gnuplot --
- Include these lines of code at the top of your c or c++ program
- to use gnuplot for graphical output. Call open_gnupipe(), then
- use gnuplot("command string") to generate plots. Finally call
- close_gnupipe() before exiting the program. This routine
- assumes that your path includes the gnuplot program.
- (Dr. Mark C. Fair)
- ***************************************************************/
- #include <stdlib.h>
- #include <stdio.h>
- #include <string.h>
- FILE *gnu_pipe;
- void open_gnupipe(){
- if (( gnu_pipe = popen("gnuplot -persist","w")) == NULL){
- perror("popen failed");
- exit(1); }
- }
- void gnuplot(char *cmdstr){
- char cmdstr1[strlen(cmdstr)+1];
- sprintf(cmdstr1,"%s\n",cmdstr);
- fputs(cmdstr1,gnu_pipe);
- fflush(gnu_pipe);
- }
- void close_gnupipe(){
- if (pclose(gnu_pipe) == -1)
- printf("Problem closing communication to Gnuplot\n");
- }
- /***************************************************************
- -- End of code to use gnuplot --
- ***************************************************************/
- #include <iostream>
- #include <fstream>
- #include <cmath>
- #include <cassert>
- using std::cout;
- using std::cin;
- using std::endl;
- using std::ofstream;
- void gravrk(double*,double,double*,double*);
- void rk4( double*,int,double,double,
- void (*derivsRK)(double*,double,double*,double*),
- double*);
- int main()
- {
- double param[5];
- param[0] = 0.5;
- param[1] = 1;
- param[2] = .05;
- param[3] = .1;
- param[4] = .01;
- double state[4];
- double nState = 4;
- double tau = .1;
- double time = 0;
- int nStep = 400;
- double *tplot,*rplot,*thplot;
- tplot = new double[nStep];
- rplot = new double[nStep];
- thplot = new double[nStep];
- double r[2],v[2];
- v[0] = 0;
- v[1] = 0;
- cout << "enter in the initial height(cm): ";
- cin >> r[0];
- cout << "enter in the initial angle(rads): ";
- cin >> r[1];
- state[0] = r[0];
- state[1] = r[1];
- state[2] = v[0];
- state[3] = v[1];
- for(int iStep = 0; iStep < nStep; ++iStep){
- rplot[iStep] = r[0]; // Record position for plotting
- thplot[iStep] = r[1];
- tplot[iStep] = time;
- rk4( state, nState, time, tau, gravrk, param );
- r[0] = state[0]; r[1] = state[1]; // 4th order Runge-Kutta
- v[0] = state[2]; v[1] = state[3];
- time += tau;
- }
- ofstream fout("12.txt");
- for(int i = 0; i < nStep; ++i)
- {
- fout << rplot[i] << " " << thplot[i] << " " << tplot[i] << endl;
- }
- open_gnupipe();
- gnuplot((char *)"set term x11 3 font 'Times,14'; \
- set title 'wilberforce pendulum' font 'Times,18'; \
- set grid; \
- set xlabel 'Time(s)' font 'Times,18'; \
- set ylabel 'z(cm) theta(radians)' font 'Times,18'; \
- plot '12.txt' using 3:1 with lines title 'Z(t)'; \
- replot '12.txt' using 3:2 with lines title 'theta(t)';\
- set term postscript eps color 'Times' 18; \
- set output '12.eps'; \
- replot");
- close_gnupipe();
- delete [] tplot,rplot,thplot;
- return 0;
- }
- void gravrk(double x[], double t, double param[], double deriv[]) {
- // Returns right-hand side of Kepler ODE; used by Runge-Kutta routines
- // Inputs
- // x State vector [r(1) r(2) v(1) v(2)]
- // t Time (not used)
- // param mass, I, k, delta, e
- // Output
- // deriv Derivatives [dr(1)/dt dr(2)/dt dv(1)/dt dv(2)/dt]
- //* Compute acceleration
- double mass = param[0];
- double I = param[1];
- double k = param[2];
- double delta = param[3];
- double e = param[4];
- double r1 = x[0], r2 = x[1]; // Unravel the vector s into
- double v1 = x[2], v2 = x[3]; // position and velocity
- double accelz = -(k/mass)*r1 - .5*(1.0/mass)*e*r2; // Gravitational acceleration
- double acceltheta = -(delta/I)*r2 - .5*(1.0/I)*e*r1;
- cout << "v1: " << v1 << endl; //<< "v2: " << v2 << endl;
- //cout << "accelz: " << accelz << endl << "acceltheta: " << acceltheta << endl;
- //* Return derivatives [dr[1]/dt dr[2]/dt dv[1]/dt dv[2]/dt]
- deriv[0] = v1; deriv[2] = v2;
- deriv[1] = accelz; deriv[3] = acceltheta;
- }
- void rk4(double x[], int nX, double t, double tau,
- void (*derivsRK)(double x[], double t, double param[], double deriv[]),
- double param[]) {
- // Runge-Kutta integrator (4th order)
- // Inputs
- // x Current value of dependent variable
- // nX Number of elements in dependent variable x
- // t Independent variable (usually time)
- // tau Step size (usually time step)
- // derivsRK Right hand side of the ODE; derivsRK is the
- // name of the function which returns dx/dt
- // Calling format derivsRK(x,t,param,dxdt).
- // param Extra parameters passed to derivsRK
- // Output
- // x New value of x after a step of size tau
- double *F1, *F2, *F3, *F4, *xtemp;
- F1 = new double [nX]; F2 = new double [nX];
- F3 = new double [nX]; F4 = new double [nX];
- xtemp = new double [nX];
- //* Evaluate F1 = f(x,t).
- (*derivsRK)( x, t, param, F1 );
- //* Evaluate F2 = f( x+tau*F1/2, t+tau/2 ).
- double half_tau = 0.5*tau;
- double t_half = t + half_tau;
- for(int i=0; i< nX; ++i)
- xtemp[i] = x[i] + half_tau*F1[i];
- (*derivsRK)( xtemp, t_half, param, F2 );
- //* Evaluate F3 = f( x+tau*F2/2, t+tau/2 ).
- for(int i=0; i< nX; ++i)
- xtemp[i] = x[i] + half_tau*F2[i];
- (*derivsRK)( xtemp, t_half, param, F3 );
- //* Evaluate F4 = f( x+tau*F3, t+tau ).
- double t_full = t + tau;
- for(int i=0; i< nX; ++i)
- xtemp[i] = x[i] + tau*F3[i];
- (*derivsRK)( xtemp, t_full, param, F4 );
- //* Return x(t+tau) computed from fourth-order R-K.
- for(int i=0; i< nX; ++i)
- x[i] += tau/6.*(F1[i] + F4[i] + 2.*(F2[i]+F3[i]));
- delete [] F1, F2, F3, F4, xtemp;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement