/*                                                                   */
/*  5/10/95   10/10/96       09/10/97                                */
/*dt=.001
itmax=300
ittmax=100 
omega0=.0
domegasdt=3

A0=1;
echy=1;
omega=1;
lbd=.150;
k=(2*Pi)*(2*Pi);
u=.001;
up=.0;
*/
#include <stdio.h>
#include <stdlib.h>
#include "math.h"
#include "rgraph.h"
#include "param.c"
#ifdef powerc
#include <sioux.h>
#endif
/*#include "unix.h"*/
double dt,t;
double* theta=NULL;
double A0,omega0,domegasdt,lbd,k,omega;
double Pi;
FILE *f,*g,*fe;
int nx,NL;
void axes(float *xbas,float *x0,int *nx,float *dx,
      float *ybas,float *xhaut,float *yhaut,int *npix,int *npiy);
double FF1(double x,double y,double z);
double FF2(double x,double y,double z);
      


void main(void)
{
int i,it,itt,itmax,ittmax;
int npix,npiy;
int nbt;
float fdt;
float x,x0,xbas,ybas,xhaut,yhaut,tmax;
float xp1,yp1,xp2,yp2,echy;
float xu1,yup1,xu2,yup2;
float us;
int icoul;
double u,up;
double du1,dup1,du2,dup2,du3,dup3,du4,dup4;
int itom,itomax;

/* Place the window¹s top left corner at (5,40).      */
 #ifdef powerc
SIOUXSettings.toppixel =80;
SIOUXSettings.leftpixel =20; 
SIOUXSettings.rows = 10;
//SIOUXSettings.columns = 40;
 #endif
 
/* ouverture des fichiers de sauvegarde  */
g = fopen("F.IN", "r");
printf("Bench E24 L14  Stokes .12 Lapl.17 \n");
printf("                                                 \n");
printf(" Resolution du pendule force                     \n");
printf("  2                                              \n");
printf(" d u        du                                   \n");
printf(" --- + lbd  -- + k sin(u) - omega sin(u) cos(u)  \n");
printf("   2        dt                                   \n");
printf(" dt                                              \n");
printf("                                                 \n");
printf("   f(u)= u     lineaire !                                \n");
/* lecture du fichier donnees */
fscanf(g,"dt=%f     \n",&fdt);
dt=fdt;
printf("pas de temps   dt=%f\n",dt);
fscanf(g,"itmax=%i          \n",&itmax);
fscanf(g,"ittmax=%i         \n",&ittmax);
printf("affichage toutes les %i iter, soit dt=%f \n",ittmax,dt*ittmax);
tmax=itmax*ittmax*dt;
printf("itmax=%d   t max=itmax*ittmax*dt= %7f  \n",itmax,tmax);
fscanf(g,"omega0=%f \n",&xp1);
omega0=xp1;
fscanf(g,"domegasdt=%f \n",&xp1);
domegasdt=xp1;
printf("omega(t)=%g+%g *t \n",omega0,domegasdt);

/* ajustement de quelques constantes */
Pi=atan(1)*4;

A0=1;
echy=.75;
omega=1;
lbd=.1;
k=(2*Pi)*(2*Pi);
//u=.0000002;
fscanf(g,"u(0)=%f \n",&xp1);
u=xp1;
//u=.02;
up=.0;



/* fin des lectures */
/* declaration */
nbt=nx+15;
theta=(double*)calloc(nbt,sizeof(double));
/* fin declaration */
/* &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*/
/*            init graphique/ trace des axes           */                  
initgraphique();
x0=0;nx=itmax*ittmax;
fdt=(float)dt;
axes(&xbas,&x0,&nx,&fdt,&ybas,&xhaut,&yhaut,&npix,&npiy); 
/* &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*/
printf(".............................................\n");
itomax=2;
for (itom=1;itom<itomax;itom++)
{
 //reffecran();

t=-dt;
//next();
for(it=0;it<itmax;it++)
{omega=omega0+domegasdt*t;
if(t<0){omega=0;}
 xp1=t;
 yp1=u/echy;
 xu1=u;
 yup1=up;
  for(itt=1;itt<ittmax;itt++)
  {
 /* Runge Kutta 4  */
   du1 =dt*FF1(u      ,up       ,t);
   dup1=dt*FF2(u      ,up       ,t);
   du2 =dt*FF1(u+du1/2,up+dup1/2,t+dt/2);
   dup2=dt*FF2(u+du1/2,up+dup1/2,t+dt/2);
   du3 =dt*FF1(u+du2/2,up+dup2/2,t+dt/2);
   dup3=dt*FF2(u+du2/2,up+dup2/2,t+dt/2);
   du4 =dt*FF1(u+du3  ,up+dup3  ,t+dt  );
   dup4=dt*FF2(u+du3  ,up+dup3  ,t+dt  );
   u =u +(du1 /6 +du2 /3 +du4 /6);
   up=up+(dup1/6 +dup2/3 +dup4/6);
   t=t+dt;
   
   
 //  u =u+ dt*FF1(u,up,t);
 //  up=up+dt*FF1(u,up,t);
  }
  printf(" %g %g %g  \n",t,u, omega*omega/k);





/* tracés */
 xp2=t;
 yp2=u/echy;
 icoul=100;
 metligne(xbas,ybas,xhaut,yhaut,npix,npiy,xp1,yp1,xp2,yp2,icoul);
 xp2=t;
 yp2=omega;
 icoul=222;
 metpoint(xbas,ybas,xhaut,yhaut,npix,npiy,xp2,yp2,icoul);
 
 xp2=t;
 us= omega/k-1 ;
 if(us<1){ 
          yp2=us/echy;
            metpoint(xbas,ybas,xhaut,yhaut,npix,npiy,xp2,yp2,icoul);}
         
 yp2=u/echy;
 
 /* trace des phases */
 xu2=u*tmax/8+tmax+6;
 xu1=xu1*tmax/8+tmax+6;
 yup2=up/echy/2/Pi;
 icoul=158;
 //metligne(xbas,ybas,xhaut,yhaut,npix,npiy,xu1,yup1,xu2,yup2,icoul);
 
 }

}
}

/*---------------------------------------------------------*/






/*---------------------------------------------------------*/
void axes(float *xbas,float *x0,int *nx,float *dx,
      float *ybas,float *xhaut,float *yhaut,int *npix,int *npiy)
{
/* trace des axes */
int i;
float xp,yp;
*xbas=*x0-1;
*ybas=-3;
*xhaut=*x0+*nx*(*dx)+.5;
*yhaut=2;
*npix= 350;
*npiy= 250;
for (i=0;i<*npix; i++)
{xp=i* *xhaut * 2 * 1.f / *npix;yp=0;
 metpoint(*xbas,*ybas,*xhaut,*yhaut,*npix,*npiy,xp,yp,i);
 xp=0;yp=i*(*dx);
 metpoint(*xbas,*ybas,*xhaut,*yhaut,*npix,*npiy,xp,yp,i);
 }
 //metligne(*xbas,*ybas,*xhaut,*yhaut,*npix,*npiy,0,0,*xhaut,0,12);getchar();
 }