/*                      
			ANFANGSWERTPROBLEM DER SCHROEDINGER GLEICHUNG
			
			Praktikum: Computational Physics, FSU Jena 
			
			Programm zur Berechnung von eindimensionalen Anfangswertproblemen der  
			Schroedinger Gleichung mit Hilfe des Crank-Nichelson Verfahrens. Die 
			Ausgabe des Betragsquadrates der Wellenfunktion erfolgt in eine Datei.  
			
			erstellt von C. Thierfelder, Jena Dezember 2002 
			
			WGS, modified on April 14, 2003

			Karsten Hannewald, modified Nov 24, 2005			
*/ 

# include <math.h> 
# include <stdio.h>
# include "cmplx.h"
# define N 1000      /*Zahl der Gitterpunkte*/

void Crank_Nichelson() 
{ 
  /* N ist die Zahl der Gitterpunkte und wird ganz oben durch define definiert*/ 
  int T=10000;              	                /* Zahl der Zeitschritte*/ 
  
  struct Complex Psi[N];	                /* Feld fuer Wellenfunktion*/ 
  struct Complex a[N],b[N];	                /* Parameter zur Wellenfunktionsberechnung*/ 
  struct Complex c,L;
  int i,j;	                                /* Zaehlvariablen*/ 
  int block=(int)(T/500);                                 /* Abstand von Zeitschritten bei Ausgabe 
						   (xfarbe erlaubt nur ca. 500*500 Pixel) */ 
  double V[N];					/* Potential*/ 
  double dx=0.001;				/* Abstand der Gitterpunkte*/ 
  double dt=0.000005;				/* Zeitschritt*/ 
  double k=30;					/* Impuls */ 
  double s=0.08;				/* Breite des Gauss  */ 
  double x0=300;				/* Anfangslage des Gauss*/ 
  double Integral_anfang=0;		        /* Flaechen unter Wellenfunktion */ 
  double Integral_kontrolle=0; 
  
  FILE *datei_ptr;				/* FILE-Pointer*/ 
  char *name="wave.dat";		        /* Name der Datei in die ausgegeben wird*/ 
  datei_ptr=fopen(name,"w");
  
  // fprintf(datei_ptr,"Zeitentw. des Betragsquadr. der Wellenfkt. \n 500 500 \n");
  
  /* Initialisierung des Potentials */ 
  for(i=0;i<N;i++)
      {
               if ( ((i-N*0.65)<N*0.05) && ((N*0.65-i)<N*0.05) )
                  V[i]=1000.0;
          else V[i]=0.0;
               //V[i]=fabs(i-N/2)*fabs(i-N/2)*0.1;
      }
  
  /* Initialisierung der Anfangswellenfunktion */ 
  for(i=1;i<N-1;i++) 
    { 
      Psi[i].re=exp(-((i*dx-x0*dx)*(i*dx-x0*dx))/(4.0*s*s))*cos(k*i*dx)/
	sqrt(sqrt(2.0*3.141)*s); 
      Psi[i].im=exp(-((i*dx-x0*dx)*(i*dx-x0*dx))/(4.0*s*s))*sin(k*i*dx)/
	sqrt(sqrt(2.0*3.141)*s); 
      Integral_anfang += Betrag(Psi[i])*Betrag(Psi[i])*dx; 
    } 
  
  Psi[N-1].re=0;   
  Psi[N-1].im=0; 
  Psi[0].re=0;   
  Psi[0].im=0;
  
  Integral_kontrolle=Integral_anfang;  
  
  /* Berechnung des Parameters a */ 
  
  a[1].re=2+dx*dx*V[1]; 
  a[1].im=-2*dx*dx/dt; 
  
  for(i=2;i<N;i++) 
    { 
      a[i].re=2+dx*dx*V[i] - Re(inv(a[i-1])); 
      a[i].im=-2*dx*dx/dt  - Im(inv(a[i-1])); 
    } 
  
  /***************************Zeitentwicklungsschleife**************************************/ 
  for(j=0;j<T;j++)		 
  {
        for(i=1;i<N-1;++i)
        {
            c.re = dx*dx*V[i]+2;
            c.im = 2*dx*dx/dt;
            opeq(&L,sub(sub( mult(c, Psi[i]), Psi[i+1]),Psi[i-1]));

            if (i==1)
                opeq(&b[1],L);
            else
                opeq(&b[i],add(L,divi(b[i-1],a[i-1])));
        }
        for (i = N-2; i>0; --i)
        {
            opeq(&Psi[i],divi(sub(Psi[i+1],b[i]),a[i]));
        }
        
     
      
      /* Berechnung & Ausgabe der normierten Flaeche 
	 zur Ueberpruefung der Unitaritaet in jedem block-ten Zeitschritt */ 
      
      if (j%block==0) { 
	Integral_kontrolle=0;
	for(i=0;i<N;i++) Integral_kontrolle += Betrag(Psi[i])*Betrag(Psi[i])*dx; 
//	printf("%d  %lf\n",j,Integral_kontrolle/Integral_anfang); 
      }
      
      /* Ausgabe in Datei formatiert für Visualisierung mit xfarbe */ 
      /* nur jeder block-te Zeitschritt wird ausgegeben */ 
      
      if (j%block==0)
      { fprintf(datei_ptr, "\n\"Time=%f\"\n", dt*j);
	for(i=1;i<N;i+=2)
        {
             fprintf(datei_ptr,"%f %f \n", .0+dx*i,Betrag(Psi[i])*Betrag(Psi[i]));
           
             // fprintf(datei_ptr,"%f \n", Betrag(Psi[i])*Betrag(Psi[i]));
        }
      }
    }
  fclose(datei_ptr); 
} 

int main() 
{	 
  printf("Iterationsschritt und normierte Flaeche unter Wellenfkt. \n"); 
  Crank_Nichelson(); 
  return 0; 
} 
/*Definition einer Struktur fuer komplexe Zahlen*/ 
struct Complex  
{ 
  double re; 
  double im; 
}; 

void Crank_Nichelson();	                            /****Loesungsroutine****/ 
struct Complex mult(struct Complex,struct Complex); /*Multiplikation komplexer Zahlen*/ 
struct Complex divi(struct Complex,struct Complex); /*Division komplexer Zahlen*/ 
struct Complex add(struct Complex,struct Complex);  /*Addition komplexer Zahlen*/ 
struct Complex sub(struct Complex,struct Complex);  /*Subtraktion komplexer Zahlen*/ 
struct Complex inv(struct Complex);	            /*Inverses einer komplexen Zahlen*/ 
double Re(struct Complex);                          /*Realteil einer komplexen Zahlen*/ 
double Im(struct Complex);                          /*Imaginaerteil einer komplexen Zahlen*/ 
double Betrag(struct Complex);                      /*Betrag einer komplexen Zahlen*/ 
void opeq(struct  Complex*, struct Complex);

void opeq(struct Complex* a, struct Complex b)
{
    (*a).re = b.re;
    (*a).im = b.im;
}

/*Komplexe Multiplikation*/ 
struct Complex mult(struct Complex a, struct Complex b) 
{ 
  struct Complex ergebnis; 
  
  ergebnis.re = a.re*b.re-a.im*b.im; 
  ergebnis.im = a.im*b.re+a.re*b.im; 
  
  return ergebnis; 
} 

/*Komplexe Division*/ 
struct Complex divi(struct Complex a, struct Complex b) 
{ 
  struct Complex ergebnis; 
  
  ergebnis.re = (a.re*b.re+a.im*b.im)/(b.re*b.re+b.im*b.im); 
  ergebnis.im = (a.im*b.re-a.re*b.im)/(b.re*b.re+b.im*b.im); 
  
  return ergebnis; 
} 

/*Komplexe Addition*/ 
struct Complex add(struct Complex a, struct Complex b) 
{ 
  struct Complex ergebnis; 
  
  ergebnis.re = a.re+b.re; 
  ergebnis.im = a.im+b.im; 
  
  return ergebnis; 
} 

/*Komplexe Subtraktion*/ 
struct Complex sub(struct Complex a, struct Complex b) 
{ 
  struct Complex ergebnis; 
  
  ergebnis.re = a.re-b.re; 
  ergebnis.im = a.im-b.im; 
  
  return ergebnis; 
} 


/*Inverses einer komplexen Zahl*/ 
struct Complex inv(struct Complex a) 
{ 
  struct Complex ergebnis; 
  
  ergebnis.re =  a.re/(a.re*a.re+a.im*a.im); 
  ergebnis.im = -a.im/(a.re*a.re+a.im*a.im); 
  
  return ergebnis; 
} 

/*Realteil einer komplexen Zahl*/ 
double Re(struct Complex a) 
{ 
  double ergebnis; 
  
  ergebnis = a.re; 
  
  return ergebnis; 
} 

/*Imaginaerteil einer komplexen Zahl*/ 
double Im(struct Complex a) 
{ 
  double ergebnis; 
  
  ergebnis = a.im; 
  
  return ergebnis; 
} 

 
/*Betrag einer komplexen Zahl*/ 
double Betrag(struct Complex a) 
{ 
  double ergebnis; 
 
  ergebnis = sqrt(a.re*a.re + a.im*a.im); 
	 
  return ergebnis; 
}