/* 25/04/95 */	
/* 28/04/95  avec calloc*/
/* ====================================================================== */	
/* resolution de l'attarcteur de Lorentz                                  */
/*                                                                        */
/*     x' = - s * x + s * y                                               */
/*     y' = - x * z + r * x - y                                           */
/*     z' =   x * y - b * z                                               */
/*                                                                        */ 
/*                                                                        */ 
/* Biblio                                                                 */ 
/*    E.A. Jackson 'Perspectives of non linear dynamics'  Cambridge 1991  */
/*         p 143                                                          */ 
#include <stdio.h>
#include <stdlib.h>
#include "math.h"
#include "rgraph.h"
#include "param.c"
#ifdef powerc
#include <sioux.h>
#endif
/*#include "unix.h"*/
 

/* ====================================================================== */
float  s,r,b;
float* F=NULL;
float t,dt,tmax;
float xpf1,xpf2,zpf1,ypf1,ypf2,zpf2,xpf3,ypf3,zpf3;
float x,x0,xbas,ybas,xhaut,yhaut;
float f1o,f2o,f3o;
int npix,npiy;
float u1,u2,u3;
int travu;
void MainEvent(void);
void axes(float *xbas,float *x0,int *nx,float *dx,
      float *ybas,float *xhaut,float *yhaut,int *npix,int *npiy); 
void MainEvent(void);
void plotxtytzt(int ittt);
void pointcarre(int ittt);
void tracesauve( );
/* ====================================================================== */

void main(void)
{
int k,Neq,jt;
long nbt;
int nt;

/*float F[nmax][Neqmax];*/
float F1[Neqmax], F2[Neqmax], F3[Neqmax], F4[Neqmax];
float Fp1[Neqmax],Fp2[Neqmax],Fp3[Neqmax],Fp4[Neqmax];

FILE *g;

/* 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
 
printf(" exemple qui marche,                      \n");
printf("      parametres dans le fichier F.IN     \n"); 
printf("  dt=.025                                 \n");
printf("  s=3;r=26.5;b=1;                         \n");
printf(" a partir de maintenant:                  \n");
printf(" PPCMac 6100 60 23s  66 20s 80 17s 600 it \n"); 
printf(" ici:                                     \n");
printf("                                          \n");
printf("     x' = - s * x + s * y                 \n");
printf("     y' = - x * z + r * x - y             \n");
printf("     z' =   x * y - b * z                 \n");


/* ouverture du fichier de lecture     */
g = fopen("F.IN", "r");
/*dt=.025;tmax=30,s=3;r=26.5;b=1;*/
fscanf(g,"dt=%f                  \n",&dt);
printf("pas de temps   dt=%f  \n",dt);
fscanf(g,"tmax=%f                  \n",&tmax);
printf("temps max    tmax=%f          \n",tmax);
printf("temps max <3000 pour 2Mo  \n" );
printf(" 133k pour 100             \n" );
printf("parametres de structure    \n");
fscanf(g,"s=%f            \n",&s);
printf("                s=%f\n",s);
fscanf(g,"r=%f          \n",&r);
printf("                r=%f\n",r);
fscanf(g,"b=%f          \n",&b);
printf("                b=%f\n",b);
printf(" Conditions initiales  \n");
fscanf(g,"x0=%f          \n",&f1o);
printf("                x0=%f\n",f1o);
fscanf(g,"y0=%f          \n",&f2o);
printf("                y0=%f\n",f2o);
fscanf(g,"z0=%f          \n",&f3o);
printf("                z0=%f\n",f3o);
fclose(g);
Neq=3;
xpf1=0;
xpf2=0;
xpf3=0;
ypf1=0;
ypf2=0;
ypf3=0;
zpf1=0;
zpf2=0;
zpf3=0;

if(r>1){
/* calcul des points fixes */
xpf1=pow((b*(r-1)),.5);
xpf2=-xpf1;
ypf1=xpf1;
ypf2=-xpf1;
zpf1=r-1;
zpf2=r-1;
/* direction perpendiculaire au plan des points fixes */
u1=ypf1*zpf2-ypf2*zpf1;
u2=zpf1*xpf2-zpf2*xpf1;
u3=xpf1*ypf2-xpf2*ypf1;
travu=(int)(u1*f1o+u2*f2o+u3*f3o);
printf("%d \n",travu);
}
printf(" Points fixes         \n ");
printf("  xpf1=%f, ypf1=%f, zpf1=%f  \n",xpf1,ypf1,zpf1);
printf("  xpf2=%f, ypf2=%f, zpf2=%f  \n",xpf2,ypf2,zpf2);
printf("  xpf3=%f, ypf3=%f, zpf3=%f  \n",xpf3,ypf3,zpf3);

/* alloc memoire*/
 nbt = 1+ 3* (long)(tmax/dt); printf("nbre de valeurs=%d\n",nbt);
 /*scanf("nbt%f          \n",&nbt);*/

 F=(float*)calloc(nbt,sizeof(float));
/* fin declaration */
/* &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*/
/*            init graphique/ trace des axes           */                  
initgraphique();
nt = 4+(int)(tmax/dt);
x0=0;
yhaut=fmax(f1o,f2o);
yhaut=fmax(yhaut,f3o);
yhaut=yhaut*1.5;
ybas=fmin(f1o,f2o);
ybas=fmin(ybas,f3o);
ybas=fmin(-yhaut,f3o);
ybas=ybas*1.25;
printf("%f \n",yhaut);
if(r>1){
//xhaut=fmax(xhaut,xpf1);
yhaut=r*1.5;
        ybas=-r*1.25;}
axes(&xbas,&x0,&nt,&dt,&ybas,&xhaut,&yhaut,&npix,&npiy); 
printf("echelles %f %f \n",yhaut,ybas);
reffecran();
printf(" raidi? \n");
getchar();
/* &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&*/
printf(".............................................\n");
  
F1[0]=0;
F1[1]=f1o;
F1[2]=f2o;
F1[3]=f3o;
t=-dt;


   for (jt=0;jt<tmax/dt;jt++)
        {
         
			for (k=1; k<=Neq; k++)
				{
					F[3*jt+k] = F1[k];                 /*  F se remplit au fur et a mesure*/
				};
			t = t +dt;
		 /*	  printf(".");	 */
/*========================================================================================*/
   /*  a partir d ici calcul standard RK4*/
			Xprime(Fp1, F1);                 
			for (k=1; k<=Neq; k++)
				{
					F2[k] = F1[k] + Fp1[k] * dt / 2;
				};
			Xprime(Fp2, F2);
			for (k=1; k<=Neq; k++)
				{
					F3[k] = F1[k] + Fp2[k] * dt / 2;
				};
			Xprime(Fp3, F3);			
			for (k=1; k<=Neq; k++)
				{
					F4[k] = F1[k] + Fp3[k] * dt;
				};
			Xprime(Fp4, F4);
			for (k=1; k<=Neq; k++)                       /* ultime accroissement de F1*/
				{
					F1[k] = F1[k] + (Fp1[k] / 6 + Fp2[k] / 3 + Fp3[k] / 3 + Fp4[k] / 6) * dt;
				};
	if(jt>0){	plotxtytzt(jt);}
		if(r>1){pointcarre(jt);}		
	 	/*  printf(" t= %7f %f %f  %f %f  \n",t,F[3*jt+1],F1[1],F[3*jt+2],F1[2]) ;  */
	 	
	}
/*========================================================================================*/	
	    printf(" FIN ");
	    getchar();
	    tracesauve();
	   ExitToShell();
	}
/*========================================================================================*/
void pointcarre(jt)
 { extern float u1,u2,u3;
 extern int travu;
  extern float x,x0,xbas,ybas,xhaut,yhaut;
  extern int npix,npiy;
  int icoul;
  float xp1,xp2,yp1,yp2;
  extern float* F;
  extern float t,dt,tmax;
  float echy;
 int trav;
   trav=(int)(u1*F[3*(jt-1)+1]+u2*F[3*(jt-1)+2]+u3*F[3*(jt-1)+3]);
   trav=trav*(u1*F[3*(jt )+1]+u2*F[3*(jt )+2]+u3*F[3*(jt )+3]);
   if(trav<0){xp1=4*tmax/3 + F[3*(jt )+1];yp1=F[3*(jt )+2];
   metpoint(xbas,ybas,xhaut,yhaut,npix,npiy,xp1,yp1,icoul);
   }
   
  }
 
/*========================================================================================*/
void plotxtytzt(jt)
 {
  extern float x,x0,xbas,ybas,xhaut,yhaut;
  extern float f1o,f2o,f3o;
  extern int npix,npiy;
  int icoul;
  float xp1,xp2,yp1,yp2,yp3;
  extern float* F;
  extern float t,dt,tmax;
  extern float s,r,b;
  float echy,lbd2;
/*trace x */
  echy=1;
  xp1=t-dt;
  yp1= F[3*(jt-1)+1]*echy;
  xp2=t;
  yp2=F[3*jt+1]*echy ;
  icoul=10; /*rouge*/
   metligne(xbas,ybas,xhaut,yhaut,npix,npiy,xp1,yp1,xp2,yp2,icoul);
/* trace y */
   echy=1;
  xp1=t-dt;
  yp1= F[3*(jt-1)+2]*echy;
  xp2=t;
  yp2=F[3*jt+2]*echy ;
  icoul=50; /* vert*/
    metligne(xbas,ybas,xhaut,yhaut,npix,npiy,xp1,yp1,xp2,yp2,icoul);
/* trace z */
   echy=1;
  xp1=t-dt;
  yp1= F[3*(jt-1)+3]*echy;
  xp2=t;
  yp2=F[3*jt+3]*echy ;
  icoul=100; /* bleu*/ 
  metligne(xbas,ybas,xhaut,yhaut,npix,npiy,xp1,yp1,xp2,yp2,icoul); 
  
/* comparaison vecteurs propres */
/* on compre a exp-bt solution evidente de dz/dt=-b z */  /* trace z */
   echy=1;
  icoul=1;
  xp1=t-dt;
  yp1=(f3o)*exp(-b*(t-dt));
  xp2=t;
  yp2=(f3o)*exp(-b*t);
  metligne(xbas,ybas,xhaut,yhaut,npix,npiy,xp1,yp1,xp2,yp2,icoul);  
 /* un autre moins trivial */
//  lbd2=(-1 - s - sqrt(-4*(1 - r)*s + pow(1 + s,2)))/2;
 // lbd2=(-1 - s + sqrt(-4*(1 - r)*s + pow(1 + s,2)))/2;
/*  xp1=t;
  yp1=-(s+lbd2)*F[3*jt+1]+r*F[3*jt+2];
  metpoint(xbas,ybas,xhaut,yhaut,npix,npiy,xp1,yp1,icoul); 
  xp1=t;
  yp3=(-(s+lbd2)*f1o+r*f2o)*exp(lbd2*(t));
  metpoint(xbas,ybas,xhaut,yhaut,npix,npiy,xp1,yp3,icoul); 
  printf("%f,%f,%f %f\n",t,yp1,yp3,yp2);*/
  
//  xp2=t;
  //yp2=F[3*jt+3];
 // metpoint(xbas,ybas,xhaut,yhaut,npix,npiy,xp2,yp2,icoul); 
  	} 
/*========================================================================================*/
void tracesauve()
 {
  extern float x,x0,xbas,ybas,xhaut,yhaut;
  extern int npix,npiy;
  FILE *f;
  int is,jt,icoul;
  float xp1,xp2,yp1,yp2;
  extern float t,dt,tmax;
  extern float* F;
 /* ouverture des fichiers de sauvegarde  */
       f = fopen("F.OUT","w");
       
		 printf(" sauvegarde  ... \n");/*scanf(" ",&is);*/
		is=0;
		for (jt=0; jt<tmax/dt-2; jt++)
		 {   
		 fprintf(f," %f %f %f %d  \n",F[3*jt+1],F[3*jt+2],F[3*jt+3],is);
		 is=2; 
		/* icoul=10; */
		//printf( " %f %f %f %f   \n",F[3*jt+1],F[3*jt+2],F[3*jt+3],t);
		/*80 100 bleu*/
        /* metligne(xbas,ybas,xhaut,yhaut,npix,npiy,xp,yp,xp2,yp2,icoul);*/
		 }
		 printf("  ... dans F.OUT  \n");
		fclose(f);
}
/*========================================================================================*/

/*---------------------------------------------------------*/
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-0.1;
//*ybas=-15;
*xhaut=*x0+*nx*(*dx);
//*yhaut=10;
*npix= 550;
*npiy= 350;
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);
 }
 }