/****** * * A C version of Sean Lynn's TAKEOFF2 FORTRAN Program * implemented by Pete MacMillin, with refinements to the numerics * done for W.H. Mason: Lynn(1993-94), MacMillin(Fall 1994) * * the input must be defined in a file: takeoff.in * the output is placed in takeoff.out * *****/ #pragma options(!assign_registers) #include #include #include /* Functions */ void Input_Data(double *rho, double *W, double *Sw, double *CLmax, double *CLgrd, double *CLair, double *CDgrd, double *CDair, double *MUgrd, double *MUbrk, double *lambda, double *k, double *t_brake, double *h_obst, int *Neng, double *T, double *V, double *dt, double *t_rot, int *print, char *case_title); void Qinterp(double*, double*, double*, double*, double*); void Take_off(double, double, double, double, double, double, double, double, double, double, double, double, double, double, int, double*, double*, double, double, double, char*, int, FILE*); void Ground_Run(double, double, double, double*, double*, double*, double*, int, FILE*); void Ground(double, double*, double*, double*); void Ground_time(double, double, double, double*, double*, double*, double*, int, FILE*); void Climb_Out(double, double, double, double, double, double*, double*, double*, double*, double*, double*, int, FILE*); void Air(double, double*, double*, double*); double Lin_interp(double, double, double, double, double); void rkqs(double*, double*, int, double*, double, double, double*, double*, double*, double*, void (*)(double, double*, double*, double*)); void rkck(double y[], double dydx[], int n, double x, double h, double yout[], double yerr[], void (*derivs)(double, double[], double[], double[]), double*); double *dvector(int nl, int nh); void free_dvector(double *v, int nl, int nh); void bomb(char*); static double at, bt, ct, maxarg1, maxarg2; #define TINY 1.0e-20; #define PYTHAG(a, b) ((at = fabs(a)) > (bt = fabs(b)) ? \ (ct = bt / at, at * sqrt(1.0 + ct*ct)) : \ (bt ? (ct = at/bt, bt * sqrt(1.0 + ct*ct)) : 0.0)) #define MAX(a, b) (maxarg1 = (a), maxarg2 = (b), \ (maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2)) #define SIGN(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) #define CA 0.0003 #define CB 1.0e-9 #define PIO2 1.57079632679490 /* Added for ODE solver, pem 7/27/93 */ #define SAFETY 0.9 #define PGROW -0.2 #define PSHRNK -0.25 #define ERRCON 1.39e-4 void main() { FILE *out; double rho, W, Sw, CLmax, CLgrd, CLair, CDgrd, CDair, MUgrd, MUbrk, lambda, k, t_brake, h_obst, T[3], V[3], dt, t_rot, g; int Neng, i, print; char case_title[80]; g = 32.174; dt = 1; Input_Data(&rho, &W, &Sw, &CLmax, &CLgrd, &CLair, &CDgrd, &CDair, &MUgrd, &MUbrk, &lambda, &k, &t_brake, &h_obst, &Neng, T, V, &dt, &t_rot, &print,case_title); out = fopen("takeoff.out","w"); Take_off(rho, W, Sw, CLmax, CLgrd, CLair, CDgrd, CDair, MUgrd, MUbrk, lambda, k, t_brake, h_obst, Neng, T, V, dt, t_rot, g, case_title, print, out); fclose(out); } void Input_Data(double *rho, double *W, double *Sw, double *CLmax, double *CLgrd, double *CLair, double *CDgrd, double *CDair, double *MUgrd, double *MUbrk, double *lambda, double *k, double *t_brake, double *h_obst, int *Neng, double *T, double *V, double *dt, double *t_rot, int *print, char *case_title) { FILE *in; if((in = fopen("takeoff.in","r")) == NULL) bomb("Failure in opening 'takeoff.in'"); fgets(case_title,80,in); fscanf(in, "%lf%*[^\n]\n%lf%*[^\n]\n%lf%*[^\n]\n",rho,W,Sw); fscanf(in, "%lf%*[^\n]\n%lf%*[^\n]\n%lf%*[^\n]\n%lf%*[^\n]\n", CLmax, CLgrd, CLair, CDgrd); fscanf(in, "%lf%*[^\n]\n%lf%*[^\n]\n%lf%*[^\n]\n%lf%*[^\n]\n", CDair, MUgrd, MUbrk, lambda); fscanf(in, "%lf%*[^\n]\n%lf%*[^\n]\n%lf%*[^\n]\n%d%*[^\n]\n", k, t_brake, h_obst, Neng); fscanf(in, "%lf %lf %lf%*[^\n]\n%lf %lf %lf%*[^\n]\n", T, T+1, T+2, V, V+1, V+2); fscanf(in, "%lf%*[^\n]\n%d%*[^\n]\n", t_rot, print); fclose(in); } void Take_off(double rho, double W, double Sw, double CLmax, double CLgrd, double CLair, double CDgrd, double CDair, double MUgrd, double MUbrk, double lambda, double k, double t_brake, double h_obst, int Neng, double *T, double *V, double dt, double t_rot, double g, char *case_title, int print, FILE *out) { double Vr, t0, t1, t2, rparm[11], t, Tr, Xr, Tlo, Xlo, Vlo, Tobs, Xobs, Uobs, Vobs, t_guess1, t_guess2, t_guess, Xcrit, Vcrit, Tcrit1, Tcrit2, Tdecs, Xdecs, Vdecs, Tbfl, Power_frac, Tnext, Vlo_oei, Xobs_oei, Vobs_oei, Tstop, Xstop, Xdif1, Xdif2, y[4]; int i, n; /* Calculate fraction of power remaining and thrust coefficients */ Power_frac = (double) (Neng-1)/Neng; Qinterp(V, T, &t0, &t1, &t2); /* Calculate rotation (take-off) velocity */ Vr = k*sqrt((2*W)/(rho*Sw*CLmax)); /* Print out stuff so far */ if (print) { fprintf(out, "\n<<<<<<<<<<<< Input Data >>>>>>>>>>>>>>\n"); fprintf(out, "%s\n", case_title); fprintf(out, "Atmospheric density = %lf (sl/ft^2)\n",rho); fprintf(out, "Aircraft Weight = %lf (lbs)\n", W); fprintf(out, "Wing area = %lf (ft^2)\n", Sw); fprintf(out, "Maximum lift coefficient = %lf\n", CLmax); fprintf(out, "Ground lift coefficient = %lf\n", CLgrd); fprintf(out, "Air lift coefficient = %lf\n", CLair); fprintf(out, "Ground drag coefficient = %lf\n", CDgrd); fprintf(out, "Air drag coefficient = %lf\n", CDair); fprintf(out, "Rolling friction coeff. = %lf\n", MUgrd); fprintf(out, "Braking friction coeff. = %lf\n", MUbrk); fprintf(out, "Angle of thrust = %lf (rad)\n", lambda); fprintf(out, "Stall margin (k) = %lf\n", k); fprintf(out, "Descision time = %lf (sec)\n", t_brake); fprintf(out, "Obstacle height = %lf (ft)\n", h_obst); fprintf(out, "OEI power remaining = %lf\n\n", Power_frac); fprintf(out, "Thrust = %lf + %lf * V + %lf * V^2\n\n\n", t0, t1, t2); } /* All engines operating take-off */ /* Ground run */ /* Set up constants */ rparm[0] = rho; rparm[1] = W; rparm[2] = Sw; rparm[3] = CLgrd; rparm[4] = CDgrd; rparm[5] = MUgrd; rparm[6] = t0; rparm[7] = t1; rparm[8] = t2; rparm[9] = lambda; rparm[10] = g; /* Print out header */ if (print==2) { fprintf(out, " Normal Takeoff >>>>>>>>>>>>>>\n"); fprintf(out, "Time (s) x-dist (ft) x-Vel (ft/s) y-dist (ft) y-Vel (ft/s)\n"); } /* Call ground run routine */ Ground_Run(0,0,0, &Tr, &Xr, &Vr, rparm, print, out); t = Tr; y[0] = Xr; y[1] = Vr; y[2] = y[3] = 0; /* Call take-off rotation routine */ Tlo = Tr + t_rot; Ground_time(Tr, Xr, Vr, &Tlo, &Xlo, &Vlo, rparm, print, out); /* Climb phase */ /* change rparms for climbing CL and CD */ rparm[3] = CLair; rparm[4] = CDair; Climb_Out(Tlo, Xlo, Vlo, 0, 0, &Tobs, &Xobs, &Uobs, &h_obst, &Vobs, rparm, print, out); Vobs = sqrt(Uobs*Uobs + Vobs*Vobs); /* Output normal takeoff data */ fprintf(out, "n<<<<<<<<<<<<<< Normal Take-off Summary >>>>>>>>>>>>>>>\n"); fprintf(out, "Rotation Velocity = %8.3lf (ft/s)\n",Vr); fprintf(out, "Lift-off velocity = %8.3lf (ft/s)\n",Vlo); fprintf(out, "Velocity over obstacle = %8.3lf (ft/s)\n",Vobs); fprintf(out, "Rotation distance = %8.3lf (ft)\n",Xr); fprintf(out, "Lift-off distance = %8.3lf (ft)\n",Xlo); fprintf(out, "Distance to obstacle = %8.3lf (ft)\n",Xobs); fprintf(out, "Rotation time = %8.3lf (s)\n",Tr); fprintf(out, "Lift-off time = %8.3lf (s)\n",Tlo); fprintf(out, "Time to obstacle = %8.3lf (s)\n\n",Tobs); /* Find Balanced Field Length */ Tcrit2 = 0.75*Tr; print = 1; /* Set up for ground roll from t=0 to t=Tcrit */ rparm[3] = CLgrd; rparm[4] = CDgrd; rparm[5] = MUgrd; rparm[6] = t0; rparm[7] = t1; rparm[8] = t2; y[0] = y[1] = 0; /* Call ground run routine */ Ground_time(0,0,0, &Tcrit2, &Xcrit, &Vcrit, rparm, print, out); /* Reduce thrust and go from Vcrit to Vr */ rparm[6] = t0*Power_frac; rparm[7] = t1*Power_frac; rparm[8] = t2*Power_frac; Ground_Run(Tcrit2, Xcrit, Vcrit, &t, y, &Vr, rparm, print, out); /* Do rotation roll */ t += t_rot; Ground_time(t - t_rot, y[0], Vr, &t, y, &Vlo_oei, rparm, print, out); /* Climb phase, OEI */ rparm[3] = CLair; rparm[4] = CDair; Climb_Out(t, y[0], Vlo_oei, 0, 0, &t, &Xobs_oei, &Uobs, &h_obst, &Vobs_oei, rparm, print, out); Vobs_oei = sqrt(Uobs*Uobs + Vobs_oei*Vobs_oei); /* Or should we have braked? */ rparm[3] = CLgrd; rparm[4] = CDgrd; /* Ground run for t_brake (Vcrit to Vdecs) */ Tdecs = Tcrit2 + t_brake; Ground_time(Tcrit2, Xcrit, Vcrit, &Tdecs, &Xdecs, &Vdecs, rparm, print, out); /* Ground run from Vdecs to V=0 */ /* NO thrust */ rparm[6] = rparm[7] = rparm[8] = 0; rparm[5] = MUbrk; y[1] = 0; Ground_Run(Tdecs, Xdecs, Vdecs, &Tstop, &Xstop, y+1, rparm, print, out); Xdif2 = Xstop - Xobs_oei; /* Loop until converged */ Tnext = 0.95*Tr; do { Xdif1 = Xdif2; Tcrit1 = Tcrit2; Tcrit2 = Tnext; /* Set up for ground roll from t=0 to t=Tcrit */ rparm[3] = CLgrd; rparm[4] = CDgrd; rparm[5] = MUgrd; rparm[6] = t0; rparm[7] = t1; rparm[8] = t2; y[0] = y[1] = print = 1; /* Call ground run routine */ Ground_time(0,0,0, &Tcrit2, &Xcrit, &Vcrit, rparm, print, out); /* Reduce thrust and go from Vcrit to Vr */ rparm[6] = t0*Power_frac; rparm[7] = t1*Power_frac; rparm[8] = t2*Power_frac; Ground_Run(Tcrit2, Xcrit, Vcrit, &t, y, &Vr, rparm, print, out); /* Do rotation roll */ t += t_rot; Ground_time(t-t_rot, y[0], Vr, &t, y, &Vlo_oei, rparm, print, out); /* Climb phase, OEI */ rparm[3] = CLair; rparm[4] = CDair; Climb_Out(t, y[0], Vlo_oei, 0, 0, &t, &Xobs_oei, &Uobs, &h_obst, &Vobs_oei, rparm, print, out); Vobs_oei = sqrt(Uobs*Uobs + Vobs_oei*Vobs_oei); /* Or should we have braked? */ rparm[3] = CLgrd; rparm[4] = CDgrd; /* Ground run for t_brake (Vcrit to Vdecs) */ Tdecs = Tcrit2 + t_brake; Ground_time(Tcrit2, Xcrit, Vcrit, &Tdecs, &Xdecs, &Vdecs, rparm, print, out); /* Ground run from Vdecs to V=0 */ /* NO thrust */ rparm[6] = rparm[7] = rparm[8] = 0; rparm[5] = MUbrk; y[1] = 0; Ground_Run(Tdecs, Xdecs, Vdecs, &Tstop, &Xstop, y+1, rparm, print, out); Xdif2 = Xstop - Xobs_oei; Tnext = Lin_interp(0, Xdif2, Xdif1, Tcrit2, Tcrit1); } while(fabs(Xdif2) > 1e-2); /* Write OEI take-off data */ fprintf(out, "n<<<<<<<<<<<<<<<<<<<<<<< OEI Take-Off Summary >>>>>>>>>>>>>>>>>>>>\n"); fprintf(out, "Critical Velocity = %8.3lf (ft/s)\n",Vcrit); fprintf(out, "Decision Velocity = %8.3lf (ft/s)\n",Vdecs); fprintf(out, "Velocity over obstacle = %8.3lf (ft/s)\n",Vobs_oei); fprintf(out, "Critical Distance = %8.3lf (ft)\n",Xcrit); fprintf(out, "Decision Distance = %8.3lf (ft)\n",Xdecs); fprintf(out, "Balanced Field Length = %8.3lf (ft)\n",Xstop); fprintf(out, "Critical Time = %8.3lf (sec)\n",Tcrit2); fprintf(out, "Decision Time = %8.3lf (sec)\n",Tdecs); fprintf(out, "OEI Takeoff Time = %8.3lf (sec)\n",t); } /***********************************************************************/ void Ground_Run(double t0, double x0, double v0, double *tf, double *xf, double *vf, double *rparm, int print, FILE *out) { double y[2], yscal[2], y_old[2], dydx[2], t, dt, dt_did, dt_next, err; int n, i; /* Set number of equations, initial x, x', y, y', and t */ n = 2; y[0] = x0; y[1] = v0; t = t0; yscal[0] = yscal[1] = 1; err = 2.5; dt = 1.0; /* Step through time from V0 to Vf */ do { if (print==2) fprintf(out, "%7.3lf %10.3lf %10.3lf %10.3lf %10.3lf\n", t, y[0], y[1], 0.0, 0.0); for (i = 0; i <= 1; i++) y_old[i] = y[i]; Ground(t, y-1, dydx-1, rparm); rkqs(y-1, dydx-1, n, &t, dt, err, yscal-1, &dt_did, &dt_next, rparm, Ground); dt = dt_next; } while ((y[1]-v0)*(y[1]-v0) < (*vf-v0)*(*vf-v0)); /* Interpolate to find time and distance at Vf */ t = *tf = Lin_interp(*vf, y_old[1], y[1], t-dt_did, t); y[0] = *xf = Lin_interp(*vf, y_old[1], y[1], y_old[0], y[0]); y[1] = *vf; if (print==2) fprintf(out, "%7.3lf %10.3lf %10.3lf %10.3lf %10.3lf\n", t, y[0], y[1], 0.0, 0.0); } /***********************************************************************/ void Ground_time(double t0, double x0, double v0, double *tf, double *xf, double *vf, double *rparm, int print, FILE *out) { double y[2], yscal[2], y_old[2], dydx[2], t, dt, dt_did, dt_next, err; int n, i; /* Set number of equations, initial x, x', y, y', and t */ n = 2; y[0] = x0; y[1] = v0; t = t0; yscal[0] = yscal[1] = 1; dt = *tf - t0; err = 1000.0; Ground(t, y-1, dydx-1, rparm); rkqs(y-1, dydx-1, n, &t, dt, err, yscal-1, &dt_did, &dt_next, rparm, Ground); *tf = t; *xf = y[0]; *vf = y[1]; } /***********************************************************************/ void Climb_Out(double t0, double x0, double u0, double y0, double v0, double *tf, double *xf, double *uf, double *yf, double *vf, double *rparm, int print, FILE* out) { double y[4], yscal[4], y_old[4], dydx[4], t, dt, dt_did, dt_next, err; int i, n=4; /* Set initial x, x', y, y', and t */ t = t0; y[0] = x0; y[1] = u0; y[2] = y0; y[3] = v0; yscal[0] = yscal[1] = yscal[2] = yscal[3] = 1; err = 2.00; dt = (int) t + 1 - t; /* Step through time from y = 0 to y = h_obst */ do { if (print==2) fprintf(out, "%7.3lf %10.3lf %10.3lf %10.3lf %10.3lf\n", t, y[0], y[1], y[2], y[3]); for (i = 0; i <= 3; i++) y_old[i] = y[i]; Air(t, y-1, dydx-1, rparm); rkqs(y-1, dydx-1, n, &t, dt, err, yscal-1, &dt_did, &dt_next, rparm, Air); dt = dt_next; } while ((y[2] - y0)*(y[2] - y0) < (*yf - y0)*(*yf - y0)); /* Interpolate to find time, distances and velocities at y=yf */ *tf = Lin_interp(*yf, y_old[2], y[2], t-dt_did, t); y[0] = *xf = Lin_interp(*yf, y_old[2], y[2], y_old[0], y[0]); y[1] = *uf = Lin_interp(*yf, y_old[2], y[2], y_old[1], y[1]); y[3] = *vf = Lin_interp(*yf, y_old[2], y[2], y_old[3], y[3]); y[2] = *yf; if (print==2) fprintf(out, "%7.3lf %10.3lf %10.3lf %10.3lf %10.3lf\n", t, y[0], y[1], y[2], y[3]); } /***********************************************************************/ void Qinterp(double *V, double *T, double *a, double *b, double *c) { double d; d = (T[1] - T[0])/(V[1] - V[0]); *c = ((T[2] - T[1])/(V[2] - V[1]) - d)/(V[2] - V[0]); *b = d - *c*(V[0] + V[1]); *a = T[0] - d*V[0] + *c*V[0]*V[1]; } /***********************************************************************/ double Lin_interp(double x_in, double x1, double x2, double y1, double y2) { /* Calculates y(x_in) given (x1,y1) and (x2,y2) using linear interpolation */ double y; y = y1 + (y2-y1)/(x2-x1)*(x_in-x1); return y; } /***********************************************************************/ void Ground(double t, double *y, double *dydx, double *rparm) { /* Note: y[1] and dydx[1] are the first elements in these arrays */ double cnst, k0, k1, k2; /* Equation is of form du/dt = g/W(K0 + K1*u + K2*u*u), where */ /* g/W is 1/mass, and K's are constants and u is the velocity */ cnst = cos(rparm[9]) + rparm[5]*sin(rparm[9]); k0 = rparm[6]*cnst - rparm[5]*rparm[1]; k1 = rparm[7]*cnst; k2 = rparm[8]*cnst + 0.5*rparm[0]*rparm[2]*(rparm[5]*rparm[3] - rparm[4]); dydx[1] = y[2]; dydx[2] = rparm[10]/rparm[1]*(k0 + k1*y[2] + k2*y[2]*y[2]); } /***********************************************************************/ void Air(double t, double *y, double *dydx, double *rparm) { /* Note: y[1] and dydx[1] are the first elements in these arrays */ double k11, k12, k13, k21, k22, k23, gamma, v2, c1, s1; /* Equations are of the form: du/dt = g/W*[k11+k12*V+k13*V*V] */ /* dv/dt = g/W*[k21+k22*V+k23*V*V] */ gamma = atan(y[4]/y[2]); v2 = y[4]*y[4] + y[2]*y[2]; c1 = cos(rparm[9] + gamma); s1 = sin(rparm[9] + gamma); k11 = rparm[6]*c1; k12 = rparm[7]*c1; k13 = rparm[8]*c1 - 0.5*rparm[0]*rparm[2]*(rparm[4]*cos(gamma) + rparm[3]*sin(gamma)); k21 = rparm[6]*s1 - rparm[1]; k22 = rparm[7]*s1; k23 = rparm[8]*s1 + 0.5*rparm[0]*rparm[2]*(rparm[3]*cos(gamma) - rparm[4]*sin(gamma)); dydx[1] = y[2]; dydx[2] = rparm[10]/rparm[1]*(k11 + k12*sqrt(v2) + k13*v2); dydx[3] = y[4]; dydx[4] = rparm[10]/rparm[1]*(k21 + k22*sqrt(v2) + k23*v2); } /***********************************************************************\ | Fifth-order Runge-Kutta with monitoring of local truncation error to | | ensure accuracy and adjust stepsize. p.719 of Numerical Recipes. | \***********************************************************************/ void rkqs(double y[], double dydx[], int n, double *x, double htry, double eps, double yscal[], double *hdid, double *hnext, double *rparm, void (*derivs)(double, double [], double [], double[])) { int i; double errmax, h, xnew, *yerr, *ytemp, hold; yerr=dvector(1,n); ytemp=dvector(1,n); /* Set step size to the initial trial value */ h = htry; for(;;) { rkck(y, dydx, n, *x, h, ytemp, yerr, derivs, rparm); errmax = 0.0; for (i=1; i<=n; i++) errmax=MAX(errmax, fabs(yerr[i]/yscal[i])); errmax /= eps; /* If truncation error too large, reduce stepsize */ if (errmax > 1.0) { hold = h; h=SAFETY*h*pow(errmax, PSHRNK); if (h < 0.1*hold) h = 0.1*hold; xnew = (*x)+h; if (xnew == *x) bomb("Stepsize underflow in rkqs"); continue; } else { if (errmax > ERRCON) *hnext=SAFETY*h*pow(errmax,PGROW); else *hnext=5.0*h; *x += (*hdid=h); for(i=1; i<=n; i++) y[i]=ytemp[i]; break; } } free_dvector(ytemp,1,n); free_dvector(yerr,1,n); } /**********************************************************************\ | Use the fifth-order Cash-Karp Runge-Kutta method, also return an | | estimate of local truncation error using the embedded fourth-order | | method. From p.719 of Numerical Recipes | \**********************************************************************/ void rkck(double y[], double dydx[], int n, double x, double h, double yout[], double yerr[], void (*derivs)(double, double[], double[], double[]), double *rparm) { int i; static double a2=0.2, a3=0.3,a4=0.6,a5=1.0,a6=0.875,b21=0.2, b31=3.0/40.0,b32=9.0/40.0,b41=0.3,b42= -0.9,b43=1.2, b51 = -11.0/54.0, b52=2.5,b53 = -70.0/27.0,b54=35.0/27.0, b61=1631.0/55296.0,b62=175.0/512.0,b63=575.0/13824.0, b64=44275.0/110592.0,b65=253.0/4096.0,c1=37.0/378.0, c3=250.0/621.0,c4=125.0/594.0,c6=512.0/1771.0, dc5= -277.0/14336.0; double dc1=c1-2825.0/27648.0,dc3=c3-18575.0/48384.0, dc4=c4-13525.0/55296.0,dc6=c6-0.25; double *ak2,*ak3,*ak4,*ak5,*ak6,*ytemp; ak2=dvector(1,n); ak3=dvector(1,n); ak4=dvector(1,n); ak5=dvector(1,n); ak6=dvector(1,n); ytemp=dvector(1,n); for(i=1; i<=n; i++) ytemp[i]=y[i]+b21*h*dydx[i]; (*derivs)(x+a2*h,ytemp,ak2, rparm); for(i=1; i<=n; i++) ytemp[i] = y[i]+h*(b31*dydx[i]+b32*ak2[i]); (*derivs)(x+a3*h,ytemp,ak3, rparm); for(i=1; i<=n; i++) ytemp[i] = y[i]+h*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]); (*derivs)(x+a4*h,ytemp,ak4, rparm); for(i=1; i<=n; i++) ytemp[i] = y[i]+h*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]); (*derivs)(x+a5*h,ytemp,ak5, rparm); for(i=1; i<=n; i++) ytemp[i]=y[i]+h*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]); (*derivs)(x+a6*h,ytemp,ak6, rparm); for(i=1; i<=n; i++) yout[i]=y[i]+h*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]); for(i=1; i<=n; i++) yerr[i]=h*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc6*ak6[i]); free_dvector(ytemp,1,n); free_dvector(ak6,1,n); free_dvector(ak5,1,n); free_dvector(ak4,1,n); free_dvector(ak3,1,n); free_dvector(ak2,1,n); } double *dvector(nl,nh) int nl,nh; { double *v; v=(double *)calloc((nh-nl+1),(unsigned) sizeof(double)); if (!v) bomb("allocation failure in dvector()"); return v-nl; } void free_dvector(double *v, int nl, int nh) { free((char*) (v+nl)); } /**********************************************************************\ | Function to shut down the program and display a message. \**********************************************************************/ void bomb(error_text) char error_text[]; { printf("\n** UH-OH! **\n\nStopping program because:\n"); printf("\n\t%s\n\n", error_text); exit(1); }