/* SAMPLE PROGRAM CALLS ONE OF THE DE SOLVE FUNCTIONS IN SLATEC 
   ddeabm(....). Slatec documentation is in slatec/src/ddeabm.f .	*/
/* call ddeabm(&(df, neq, t, x, tout, info, rtol, atol, idid, rwork,
  				 		       lrw, rpar)) 
 df(time, X[], V[], dparams[],iparams[]) 
    is the function representing the 1st order diff eqn
 X[], and V[] are the positions and velocities at this time.		*/



#include <math.h>
#include <stdio.h>

#define DEBUG 0
#define neqn 2
#define firstCall 	0
#define scalarTol 	1
#define solveOnlyTout 	2
#define canSurpassTout 	3

#define ATOL .0005
#define RTOL .001

#define PTS 500
#define W0 1.


typedef void (*dedef)(double*,double*,double*,double*,int*);
/*typedef void (defun)(double* double*, double* double* int*); */
void df(double* t, double* X, double* V, double* dpar, int* ipar)  {
double om, al, f,x,v;
 om=dpar[0];
 al=dpar[1];
 f=dpar[2];
 x=X[0];
 v=X[1];
 V[0]=X[1];   /* X[0] is position  X[1] is velocity.  */
 V[1]=-W0*W0*sin(x)-al*v+cos(om*(*t));
} /* damped driven pendulum */



main() {
 FILE *fp, *fq, *fo;
double par[3];
double omegao, alpha, fdriv, xo, vo,t, tout;
float f;
double d,dt;
int ipar[1];
int neq, lrw,liw,status,i,j,k;
int info[15];
double x[2];
double rtol,atol;
double rwork[172]; /* 130+21*neqn are required  */
int iwork[100];
dedef F;
F=&df;
 fp=fopen("p.dat","w");	/* open 3 data files to record:		*/
 fq=fopen("q.dat","w");	/* 	x(t), v(t), v(x).		*/ 
 fo=fopen("o.dat","w");
rtol=RTOL;atol=ATOL;
neq=neqn; lrw=172; liw=100;

puts("omega_o, alpha (dampening), xo (init pos/amplitude), vo, \n and f (force val; if applicable; set f=0 to automatic f=fstable)  ?\n");
scanf("%f", &f); omegao=1.*f;
scanf("%f", &f); alpha=1.*f;
scanf("%f", &f); xo=1.*f;
scanf("%f", &f); vo=1.*f;
scanf("%f", &f); fdriv=1.*f;
if(fdriv==0.) {
 /* find effective amplitude */
 d=.5*omegao*omegao*xo*xo+.5*vo*vo;
 d=sqrt(2.*d/omegao/omegao);
 fdriv=omegao*omegao*d*alpha/2./3.14159;
};
printf("Fdriv:=%g \n", fdriv);
par[0]=omegao;
par[1]=alpha;
par[2]=fdriv;

x[0]=xo;
x[1]=vo;

puts("Final time? \n");
scanf("%f", &f); tout=1.*f;
t=0.;
dt=tout/PTS;
printf("Time interval is:(%g,%g)\n", t,tout);


info[firstCall]=0;
info[scalarTol]=0;
info[solveOnlyTout]=1;
info[canSurpassTout]=0;


ddeabm(F,&neq,&t,x,&tout,info,&rtol,&atol,&status,rwork,&lrw,iwork,&liw,par,ipar);


for(i=1,j=-1,k=0;(i<=1);j=k) 
{ if (DEBUG) {printf(" S<%d",status);puts("XOX");}
  k=(int) (t/dt);
  if ((k-j)>0) {
  	fprintf(fo,"%g     %g \n", x[0], x[1]);
 	fprintf(fq,"%g     %g \n", t, x[0]);
 	fprintf(fp,"%g     %g \n", t, x[1]);	
 }
 
  i=status;
  if (i==-2) info[firstCall]=1;
   else if(i<0) {puts("Error happened. exiting.\n");break;}
  if (i<=1)  
ddeabm(F,&neq,&t,x,&tout,info,&rtol,&atol,&status,rwork,&lrw,iwork,&liw,par,ipar);
 
}

fclose(fp);
fclose(fq);
fclose(fo);
}