/* author: V. P. Singh */


/*  Time evolution of wavefunction of cloud 1 and 2 for each sample     */



#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "functions.c"


const int xlen = 200 ;
const int ylen = 200 ;


const double J0 =  1.0 ;
const double  U =  0.156 ;



double freq_x_Hz = 11;
double freq_y_Hz = 11;

double trap_scale = 4.6208e-6;          /*0.5 m*(2Pi l )^2/J  */

double J12 = J0;

const double  eps=1.e-8 ;                /* precision for Runge-Kutta ODE integration */

/* Prefactor turn off parameters  */

double ti = 0.0 ;
double tf;



double J12_off(double t)
{
  return J12*exp(-1.*t/tf );
}




double gdamp;



/* Copying real and imag parts of the derivative to "dxdt" array owing to GPE equation of motion for each lattice site*/
void dervis(double t, double *xfield, double *dxdt)
{
  for(int y_site=0; y_site<ylen; y_site++){
    for(int x_site=0; x_site<xlen; x_site++){


      double left_s1_re = xfield[i2d_left_re(x_site, y_site, xlen, ylen)];
      double left_s1_im = xfield[i2d_left_im(x_site, y_site, xlen, ylen)];

      double right_s1_re = xfield[i2d_right_re(x_site, y_site, xlen, ylen)];
      double right_s1_im = xfield[i2d_right_im(x_site, y_site, xlen, ylen)];

      double down_s1_re = xfield[i2d_down_re(x_site, y_site, xlen, ylen)];
      double down_s1_im = xfield[i2d_down_im(x_site, y_site, xlen, ylen)];

      double up_s1_re = xfield[i2d_up_re(x_site, y_site, xlen, ylen)];
      double up_s1_im = xfield[i2d_up_im(x_site, y_site, xlen, ylen)];


      double field_s1_re = xfield[i2d_re(x_site, y_site, xlen)];
      double field_s1_im = xfield[i2d_im(x_site, y_site, xlen)];

      double den_s1 = pow( field_s1_re, 2) + pow( field_s1_im, 2);

      


      double left_s2_re = xfield[i2d_s2_left_re(x_site, y_site, xlen, ylen)];
      double left_s2_im = xfield[i2d_s2_left_im(x_site, y_site, xlen, ylen)];

      double right_s2_re = xfield[i2d_s2_right_re(x_site, y_site, xlen, ylen)];
      double right_s2_im = xfield[i2d_s2_right_im(x_site, y_site, xlen, ylen)];

      double down_s2_re = xfield[i2d_s2_down_re(x_site, y_site, xlen, ylen)];
      double down_s2_im = xfield[i2d_s2_down_im(x_site, y_site, xlen, ylen)];

      double up_s2_re = xfield[i2d_s2_up_re(x_site, y_site, xlen, ylen)];
      double up_s2_im = xfield[i2d_s2_up_im(x_site, y_site, xlen, ylen)];


      double field_s2_re = xfield[i2d_s2_re(x_site, y_site, xlen, ylen)];
      double field_s2_im = xfield[i2d_s2_im(x_site, y_site, xlen, ylen)];

      double den_s2 = pow( field_s2_re, 2) + pow( field_s2_im, 2);



      if(  x_site == 0  )      { left_s1_re =0.0;  left_s1_im =0.0;     left_s2_re =0.0;  left_s2_im =0.0;  }
  
      if( x_site  == xlen -1 ) { right_s1_re =0.0; right_s1_im =0.0;   right_s2_re =0.0; right_s2_im =0.0;  }

      if( y_site  == 0  )       { down_s1_re=0.0;   down_s1_im=0.0;     down_s2_re=0.0;   down_s2_im=0.0;   }
  
      if( y_site  == ylen -1 ) { up_s1_re = 0.0;   up_s1_im = 0.0;     up_s2_re = 0.0;   up_s2_im = 0.0;  }



   
      
      double trap_local= trap_2d(x_site, y_site, xlen, ylen, freq_x_Hz, freq_y_Hz, trap_scale);

      dxdt[i2d_re(x_site, y_site, xlen)] = -J0*( right_s1_im + left_s1_im + up_s1_im + down_s1_im )  +    ( U*den_s1 + trap_local)*field_s1_im  - 0.5*gdamp*field_s1_re - J12_off(t)*field_s2_im ;   
      dxdt[i2d_im(x_site, y_site, xlen)] =  J0*( right_s1_re + left_s1_re + up_s1_re + down_s1_re )  -    ( U*den_s1 + trap_local)*field_s1_re  - 0.5*gdamp*field_s1_im + J12_off(t)*field_s2_re ;
      
      dxdt[i2d_s2_re(x_site, y_site, xlen, ylen)] = -J0*( right_s2_im + left_s2_im + up_s2_im + down_s2_im )  +    ( U*den_s2 + trap_local)*field_s2_im - 0.5*gdamp*field_s2_re - J12_off(t)*field_s1_im ;
      dxdt[i2d_s2_im(x_site, y_site, xlen, ylen)] =  J0*( right_s2_re + left_s2_re + up_s2_re + down_s2_re )  -    ( U*den_s2 + trap_local)*field_s2_re - 0.5*gdamp*field_s2_im + J12_off(t)*field_s1_re ;  

      
    }
  }

}







void piover2pulse(double *xfield)
{

  for(int y_site=0; y_site<ylen; y_site++){
    for(int x_site=0; x_site<xlen; x_site++){


      double field_1_re =   ( xfield[i2d_re(x_site, y_site, xlen)] + xfield[i2d_s2_re(x_site, y_site, xlen, ylen)] )/sqrt(2.0) ;
      double field_1_im =   ( xfield[i2d_im(x_site, y_site, xlen)] + xfield[i2d_s2_im(x_site, y_site, xlen, ylen)] )/sqrt(2.0) ;


      double field_2_re =   ( xfield[i2d_re(x_site, y_site, xlen)] - xfield[i2d_s2_re(x_site, y_site, xlen, ylen)] )/sqrt(2.0) ;
      double field_2_im =   ( xfield[i2d_im(x_site, y_site, xlen)] - xfield[i2d_s2_im(x_site, y_site, xlen, ylen)] )/sqrt(2.0) ;

      
      xfield[i2d_re(x_site, y_site, xlen)] = field_1_re;
      xfield[i2d_im(x_site, y_site, xlen)] = field_1_im;

   
      xfield[i2d_s2_re(x_site, y_site, xlen, ylen)] = field_2_re;
      xfield[i2d_s2_im(x_site, y_site, xlen, ylen)] = field_2_im;    

    }
  }


}




int main(int argc, char **argv)
{

  const    int Ipn  = 1; 
  const double Tmax = 2000;      /* total time for ODE integration */
  const int   Steps = 100;     /* Time Steps */
  

  int sample_number = 0;         /* Sample number of the ensemble */

  double texp = 50 ;    /* time in ms */
  
  double tunit= 0.68424;
    
  tf=  texp/tunit;

  
  
  double yrangeHalf= 5.0;  /* Corresponding to a strip of 5mum in the y direction */

  
  const double tpulse = 0.0 ;       /* quench time */
  
 
  int number_trajectories = Ipn ;

  srand(time(NULL)+sample_number);           /* initialize random seed */


  double *xfield ;
  xfield = my_vector(4*xlen*ylen);


  
  FILE *outfield, *outenergy, *outdensity, *outgfsum ;
  char filename[500];

  
  
  /* ODE Runge-Kutta (RKCB) integration */

  
  double htry, hnext, timestep, hbid ;

  
  double *dxdt, *yscal;
  dxdt  = my_vector(4*xlen*ylen);               /* derivatives, needed for the runge-kutta process.  */
  yscal = my_vector(4*xlen*ylen);               /* initialize to 1, see below. can be used to tweak the RK process, by putting different weights in the error estimate.  */
		

  
  for(int i=0; i<4*xlen*ylen; i++) yscal[i]=1;            /* initialize to 1. */

  
  
  timestep=Tmax/Steps;                                     /* Tmax is the full time. Steps the number of steps. hence... */

  htry=timestep;                                           /* attempted step size */

  hnext=htry;





  
  
  double *phase, *den, *vortex ;
    
  phase     = my_vector(xlen*Steps);
  den       = my_vector(xlen*Steps);

  vortex    = my_vector(xlen*Steps);


 
  
  double phase1, phase2, phase3, phase4, phase12, phase23, phase34, phase41, phase_site;
    


  for(int n=1; n<=Ipn; n++)      /* Ipn = total number of iterations per node */
    {
      /* reading field data from the monte-carlo initialization in harmonic trap */
      sprintf(filename, "mc_field/mc_state_field_2d_%d", ((n-1)+sample_number*Ipn));

      outfield = fopen(filename, "rb");
      if(outfield==NULL) perror("error in opening file");
      else
	{
	  for(int i=0; i<2*xlen*ylen; i++){	
	    fscanf(outfield, "%le\n", &xfield[i]);
	  }
	}
      fclose(outfield);


       /* initialization with vacuum field for second field */
      for(int i=0; i< xlen*ylen; i++){

	xfield[2*i + 2*xlen*ylen]     = gasdev()*0.5 ;
	xfield[2*i + 2*xlen*ylen + 1] = gasdev()*0.5 ;	

      }

       
      int tpulse1=0;
     
      double t=0.0;
     
      /* time integration loop */
      for(int j=0;    j<Steps;     j++){
    
	htry=timestep;
	hnext=htry;

	if(t >= tpulse && tpulse1==0 ){

	  piover2pulse(xfield);
	  tpulse1=1;
	}
	
	while(t < j*timestep){
	  if(hnext>(j*timestep-t)){htry=(j*timestep -t);} else{htry=hnext;}    
	  dervis(t, xfield, dxdt);					
	  rkqs(xfield, dxdt, 4*xlen*ylen-1, &t, htry, eps, yscal, &hbid, &hnext, dervis);
	}

	
	/* observables calculation   */

	for(int x_site=0; x_site<xlen; x_site++){

	  int y_site = ylen/2;
	
	  double field_s1_re = xfield[i2d_re(x_site, y_site, xlen)];
	  double field_s1_im = xfield[i2d_im(x_site, y_site, xlen)];
	  double field_s2_re = xfield[i2d_s2_re(x_site, y_site, xlen, ylen)];
	  double field_s2_im = xfield[i2d_s2_im(x_site, y_site, xlen, ylen)];

	  double phase_1 = atan2(field_s1_im,   field_s1_re);
	  double phase_2 = atan2(field_s2_im,   field_s2_re);
	  
	  double phase12 = phase_1  -  phase_2 ;
	  
	  if(phase12 > M_PI) phase12 -= 2*M_PI;
	  else if(phase12 < -M_PI) phase12 += 2*M_PI;
	  
	  phase[i2d(x_site, j, xlen)] = phase12;
	  
	  den[i2d(x_site, j, xlen)] =   field_s1_re*field_s1_re   +    field_s1_im*field_s1_im ;	   
   
	} 
	
      } 
    }  

  double norm_fact =1.*number_trajectories;
  double Jovhbar= 1461.48;


  sprintf(filename, "phase_exp_final/phase12_den_time_%d_%f_%d_%f.dat", sample_number, Tmax, Steps, texp);
  outgfsum = fopen(filename, "w");  
  for(int i=0; i<xlen*Steps; i++){
    fprintf(outgfsum, "%le           %le\n",    phase[i],      den[i] );
  }
  fclose(outgfsum);

 
  free(xfield);
  free(dxdt);
  free(yscal);

  free(phase);
  free(den);
  free(vortex);


  
  return 0;

 
}