/* 
   Author copyright: V. P. Singh
   Time evolution code for the Monte Carlo state with barrier parameters same as experiments. This gives the imbalance.  

*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "functions.c"
const int xlen = 508 ;
const int ylen = 33 ;
const int zlen = 33 ; 
const double J0 = 1.00 ;
const double  U = 0.89;
double freq_x_Hz = 1.6;
double freq_y_Hz = 252;
double freq_z_Hz = 250;
double trap_scale = 9.48074e-10;         
const double  eps=1.e-7 ;     
int xmid = xlen/2;
int ymid = ylen/2;
int zmid = zlen/2;
int  Lbx = xlen-3;
int  Lby = ylen-3;
double tunit = 0.0000615816;
const double ti       = 0.0 ;
const double tf       = 180.0 ;
const double Vs_finalTime = 220.0 ;
double ramp_up(double t)
{
  if(t < ti){
    return 0.0 ;
  } 
  else if( (t >= ti) && (t < tf) ){ 
    return 1.0*(t - ti)/(tf - ti) ; 
  }
  else 
    return 1.0;
}
double xx0;
double Vs, gaus_width;
/* Returnsa Gaussian potential with a half width */
double gaussian_stirrer( int x_site )
{
  double x0center  = 1.*xx0 - x_site;
  return 1.0*exp( -2.*x0center*x0center/(gaus_width*gaus_width ) );
}
/* Returns a Gaussian potential with a half width */
double gaussian_mobile( int x_site, double t, double v0, double x1, double freq){
  double om_freq =2.*M_PI*freq;   
  double x_time = v0*(t- Vs_finalTime)  +  x1*sin( om_freq*(t - Vs_finalTime) );
  double modx0 = 1.*xx0 - x_site + x_time;   
  return 1.0*exp( -2.*modx0*modx0/( gaus_width*gaus_width )  );
}
double excitation_protocol(int x_site, double t,  double v0, double v1, double freq){
  if(t < Vs_finalTime) {  return ramp_up(t)*Vs*gaussian_stirrer(x_site);   }  
  else return  Vs*gaussian_mobile(x_site, t, v0, v1, freq);
}
double *trap;
double v0, v1, simfreq;
void dervis(double t, double *xfield, double *dxdt)
{
  for(int z=0; z<zlen; z++){
    for(int y=0; y<ylen; y++){
      for(int x=0; x<xlen; x++){		
	double re = xfield[i3d_re(x, y, z, xlen, ylen)];
	double im = xfield[i3d_im(x, y, z, xlen, ylen)];
	double den = re*re  + im*im;
	double right_re = xfield[i3d_right_re(x, y, z, xlen, ylen, zlen)];
	double right_im = xfield[i3d_right_im(x, y, z, xlen, ylen, zlen)];
	double left_re = xfield[i3d_left_re(x, y, z, xlen, ylen, zlen)];
	double left_im = xfield[i3d_left_im(x, y, z, xlen, ylen, zlen)];
	double up_re = xfield[i3d_up_re(x, y, z, xlen, ylen, zlen)];
	double up_im = xfield[i3d_up_im(x, y, z, xlen, ylen, zlen)];
	double down_re = xfield[i3d_down_re(x, y, z, xlen, ylen, zlen)];
	double down_im = xfield[i3d_down_im(x, y, z, xlen, ylen, zlen)];
	double top_re = xfield[i3d_top_re(x, y, z, xlen, ylen, zlen)];
	double top_im = xfield[i3d_top_im(x, y, z, xlen, ylen, zlen)];
	double bottom_re = xfield[i3d_bottom_re(x, y, z, xlen, ylen, zlen)];
	double bottom_im = xfield[i3d_bottom_im(x, y, z, xlen, ylen, zlen)];
	if( x == 0  )      { left_re =0.0;   left_im =0.0;  } 
	if(x == xlen -1 ) { right_re =0.0;  right_im =0.0;  }
	if( y == 0  )     { down_re=0.0;   down_im=0.0;   }  
	if(y == ylen -1 ) { up_re = 0.0;   up_im = 0.0;  }
	if( z == 0  )     { bottom_re=0.0;   bottom_im=0.0;   } 
	if(z == zlen -1 ) { top_re = 0.0;   top_im = 0.0;  }
	double trap_val =trap[i3d(x, y, z, xlen, ylen)];
	double ex_pot = excitation_protocol(x, t, v0, v1, simfreq);   
	dxdt[i3d_re(x, y, z, xlen, ylen)] = -J0*(right_im + left_im + up_im + down_im + top_im + bottom_im) + U*den*im  + (trap_val  + ex_pot)*im ;
	dxdt[i3d_im(x, y, z, xlen, ylen)] = J0*(right_re + left_re + up_re + down_re + top_re + bottom_re)  - U*den*re -  (trap_val + ex_pot)*re ; 
	
      }
    }
  }
}
 

int main(int argc, char **argv)
{
  int sample_number  = atoi(argv[1]);   
  double meanmu =0.71 ;
  Vs   = atof(argv[2])*meanmu;      
  gaus_width   = atof(argv[3]);      
  v0 =  atof(argv[4]);
  v1 =  atof(argv[5]);
  double freqHz =  atof(argv[6]);
  int ncycle =  atoi(argv[7]);
  xx0  =  atof(argv[8]);
  simfreq = freqHz/(1./tunit);
  double tmove = 1.*ncycle/(freqHz*tunit);
  double Tmax = Vs_finalTime + tmove;
  int TprobeSteps= ( (int)  Vs_finalTime*1/2.) ;
  int Steps = ( (int) Tmax*1/2.) ;    
  int FinalSteps =  Steps - TprobeSteps;
  srand(time(NULL) );   
  double *xfield ;
  xfield = my_vector(2*xlen*ylen*zlen);
  FILE *outfield, *outenergy, *outdensity,  *outdensity1 ;
  char filename[500];
  trap = get_trap_3d(xlen, ylen, zlen, freq_x_Hz, freq_y_Hz,freq_z_Hz, trap_scale);
  double htry, hnext, timestep, hbid ;
  double *dxdt, *yscal;
  dxdt  = my_vector(2*xlen*ylen*zlen);                     
  yscal = my_vector(2*xlen*ylen*zlen);                     
							      
  for(int i=0; i<2*xlen*ylen*zlen; i++) yscal[i]=1;          
  timestep=Tmax/Steps;                                     
  htry=timestep;                                        
  hnext=htry;
  double   *den_Left,  *den_Right, *nv_time;
  den_Left   = my_vector(FinalSteps);
  den_Right  = my_vector(FinalSteps);
  nv_time    = my_vector(FinalSteps);
  for(int j=0; j< FinalSteps; j++){ den_Left[j]=0.0 ;   den_Right[j]=0.0 ;  }

  
  snprintf(filename,  sizeof(filename), "mc_field_htrap_%d", sample_number );
  outfield = fopen(filename, "rb");
  if(outfield==NULL) perror("error in opening file");
  else
    {
      for(int i=0; i<2*xlen*ylen*zlen; i++){	
	fscanf(outfield, "%le\n", &xfield[i]);
      }
    }
  fclose(outfield);     
  double t=0.0;
  for(int j=0; j<Steps; j++){
    htry=timestep;
    while(t < j*timestep){ 
      if(hnext>(j*timestep-t)){htry=(j*timestep -t);} else{htry=hnext;}    
      dervis(t, xfield, dxdt);					
      rkqs(xfield, dxdt, 2*xlen*ylen*zlen-1, &t, htry, eps, yscal, &hbid, &hnext, dervis);
    }
    if(j >= TprobeSteps ){
      double n_Left =0.0;
      double n_Right=0.0;
      int vortex_sum =0;
      double om_freq = 2.*M_PI*simfreq;
      double vel_time = v0*(t- Vs_finalTime)  + v1*sin( om_freq*(t - Vs_finalTime) );
      double modx0 = 1.*xx0  + vel_time;
      for(int z=0; z<zlen; z++){
	for(int y=0; y<ylen; y++){
	  for(int x=0; x<xlen; x++){
	    double re = xfield[i3d_re(x, y, z, xlen, ylen)];
	    double im = xfield[i3d_im(x, y, z, xlen, ylen)];	   
	    double density_local =  re*re  +  im*im;		
	    x < modx0 ? n_Left += density_local :  n_Right += density_local ;
	  } 		    	    
	}
      }
      den_Left[ j-TprobeSteps]  += n_Left;
      den_Right[j-TprobeSteps]  += n_Right; 
    }		
  } 
  snprintf(filename,  sizeof(filename), "den_LR_%d_%f_%f_%f_%f_%f_%f_%d.dat", sample_number, Vs/meanmu, gaus_width, xx0, v0, v1, freqHz, ncycle);
  outdensity = fopen(filename, "w");
  for(int i=0; i< FinalSteps; i++){
    fprintf(outdensity, "%le           %le\n",    den_Left[i],    den_Right[i] );
  }
  fclose(outdensity);
  free(xfield);
  free(dxdt);
  free(yscal);
  free(trap);
  free(den_Left);
  free(den_Right);
  free(nv_time);
  return 0;


}