/* 
   Author copyright: V. P. Singh
   Monte Carlo initialization for the initial cloud at a temperature T in units of J 

*/


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "functions.c"
#include "random.h"
const int xlen = 508 ;
const int ylen = 33 ;
const int zlen = 33 ; 
const double J0 = 1.00 ;
const double  U = 0.89;
const double Mu = -5.3 ;
const double  T =  0.5 ;  
const  int mc_total_steps_site = 280000;
const  int total_samples = 50 ;
const  int new_mc_updates= 1000 ;
const double sigma = 0.11 ;
double freq_x_Hz = 1.6;
double freq_y_Hz = 252;
double freq_z_Hz = 250;
double trap_scale = 9.48074e-10; 
int xmid = xlen/2;
int ymid = ylen/2;
int zmid = zlen/2;
double Ein, Eup ;
double calculate_local_energy(double *xfield, double *trap, int x, int y, int z)
{
  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 local_energy = -2.*J0*( re*(right_re + left_re + up_re + down_re + top_re + bottom_re)  + im*(right_im + left_im + up_im + down_im + top_im + bottom_im) )  + 0.5*U*den*den  +  (trap_val - Mu)*den; 
  return local_energy;
}
void mc_updated_state(double *xfield, double *trap, int mc_total_steps_per_site)
{
  for(int mc_step=0; mc_step<mc_total_steps_per_site; mc_step++)
    {
      int rand_num, new_yid, new_xid, new_zid ;  
      for(int i=0; i<xlen*ylen*zlen; i++){
	int x = randint(xlen);
	int y = randint(ylen);
	int z = randint(zlen);	  
	Ein = calculate_local_energy(xfield, trap, x, y, z);
	double gas_delta_re, gas_delta_im ;
	gas_delta_re = gasdev()*sigma;
	gas_delta_im = gasdev()*sigma;
	xfield[i3d_re(x, y, z, xlen, ylen)] += gas_delta_re;   
	xfield[i3d_im(x, y, z, xlen, ylen)] += gas_delta_im;
	Eup = calculate_local_energy(xfield, trap, x, y, z);
	if(Eup>Ein) {
	  double boltzmann_weight, random_weight;
	  random_weight = ((double)random())/(RAND_MAX + 1.0);  
	  boltzmann_weight = exp(-(Eup-Ein)/T);
	  if(random_weight > boltzmann_weight){
	    xfield[i3d_re(x, y, z, xlen, ylen)] -= gas_delta_re;  
	    xfield[i3d_im(x, y, z, xlen, ylen)] -= gas_delta_im; 
	  }
	}
      } 
    }
}


int main()
{
  
  srand(time(NULL)); 
  double *xfield;
  xfield = my_vector(2*xlen*ylen*zlen);
  double* trap;
  trap = get_trap_3d(xlen, ylen, zlen, freq_x_Hz, freq_y_Hz,freq_z_Hz, trap_scale);
  mc_updated_state(xfield, trap, mc_total_steps_site);
  FILE *outsample, *outdenxy, *outdenxz;
  char filename[300];
  for(int number_sample=0; number_sample<total_samples; number_sample++)
    {
      mc_updated_state(xfield, trap, new_mc_updates); 
      snprintf(filename,  sizeof(filename), "mc_field_htrap_%d", number_sample);
      outsample = fopen(filename, "w");
      for(int i=0; i<2*xlen*ylen*zlen; i++){
	  fprintf(outsample,"%le\n", xfield[i]);
	} 
      fclose(outsample);
    }
  free(xfield);  
  return 0;

}