/* Monte Carlo initialization for the initial cloud at a temperature T in units of J */


/* author: V. P. Singh */




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



const int x_size = 200 ;
const int y_size = 200 ;

const double J0 =  1.0 ;
const double  U =  0.144 ;



double freq_x_Hz = 13.5;
double freq_y_Hz = 13.5;

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



const double Mu = -1.54 ;
const double  T =  8.0 ;



const double sigma = 0.48 ;




const int mc_total_steps_site = 500000;           /* initial total steps per site for updating field values on each site */


const int       total_samples = 500 ;           /* total number of samples for ensemble average */
const int      new_mc_updates = 1000 ;           /* new monte carlo updates for each site for each different sample */




double initial_mean_density = fabs(Mu)/U; 


double Ein, Eup ;






/* Returns local energy for each lattice site */
double calculate_local_energy(double *xfield, int new_xid, int new_yid)
{

  double left_re = xfield[i2d_left_re(new_xid, new_yid, x_size, y_size)];
  double left_im = xfield[i2d_left_im(new_xid, new_yid, x_size, y_size)];

  double right_re = xfield[i2d_right_re(new_xid, new_yid, x_size, y_size)];
  double right_im = xfield[i2d_right_im(new_xid, new_yid, x_size, y_size)];

  double down_re = xfield[i2d_down_re(new_xid, new_yid, x_size, y_size)];
  double down_im = xfield[i2d_down_im(new_xid, new_yid, x_size, y_size)];

  double up_re = xfield[i2d_up_re(new_xid, new_yid, x_size, y_size)];
  double up_im = xfield[i2d_up_im(new_xid, new_yid, x_size, y_size)];

  
  
  if(    new_xid == 0  )   { left_re =0.0;  left_im =0.0; }
  
  if(new_xid == x_size -1 ) { right_re =0.0; right_im =0.0;  }

  if( new_yid == 0  )       { down_re=0.0;   down_im=0.0;  }
  
  if(new_yid == y_size -1 ) { up_re = 0.0;   up_im = 0.0; }
  
  

  double field_re = xfield[i2d_re(new_xid, new_yid, x_size)];
  double field_im = xfield[i2d_im(new_xid, new_yid, x_size)];

  double den = pow(field_re, 2) + pow(field_im, 2);
  
  
  return -2.*J0*( field_re*(left_re + right_re + up_re + down_re)   +field_im *(left_im + right_im + up_im + down_im) )
    
    + (U/2.)*pow( den, 2)  +
    
    ( trap_2d(new_xid, new_yid, x_size, y_size, freq_x_Hz, freq_y_Hz, trap_scale) - Mu )*den;
  
}






void mc_updated_state(double *xfield, int mc_total_steps_per_site)
{ 
 
  unsigned long long int accepted_steps=0, rejected_steps=0;

  for(int mc_step=0; mc_step<mc_total_steps_per_site; mc_step++)
    {

      int rand_num, new_yid, new_xid ;
      
      for(int i=0; i< x_size*y_size; i++){
 
	rand_num = random() % (x_size*y_size);       /* generates integer random number in between 0 and (x_size*y_size-1) */
	new_yid  = rand_num/x_size;
	new_xid  = rand_num - new_yid*x_size;
	  
	/* printf("%d %d %d\n", rand_num, new_xid, new_yid); */

	Ein = calculate_local_energy(xfield, new_xid, new_yid);

	double gas_delta_re, gas_delta_im ;

	/*double sigma_local = sigma*exp(-1.0*harmonic_potential2d(new_xid, new_yid, x_size, y_size, htrap_energyX, htrap_energyY)/(2.*T)) ; */

	gas_delta_re = gasdev()*sigma ;
	gas_delta_im = gasdev()*sigma ;
      
	xfield[i2d_re(new_xid, new_yid, x_size)] += gas_delta_re;   /* real part */
	xfield[i2d_im(new_xid, new_yid, x_size)] += gas_delta_im;   /* imag part */
      
	Eup = calculate_local_energy(xfield, new_xid, new_yid);
	/* printf("%f %f\n", Ein, Eup); */
      
	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[i2d_re(new_xid, new_yid, x_size)] -= gas_delta_re;  
		xfield[i2d_im(new_xid, new_yid, x_size)] -= gas_delta_im; 
		rejected_steps++;
		/* printf("%s\n","higher energy rejected"); */
	      }
	    else {accepted_steps++; /* printf("%s\n","higher energy accepted"); */ }
	  }
	else {accepted_steps++;  /* printf("%s\n","lower energy accepted"); */ }
    
      }

    }

  //printf("accepted steps = %llu, rejected steps =%llu, total steps=%llu, acceptance rate =%f\n", accepted_steps, rejected_steps, accepted_steps+rejected_steps, 
  	// ((double)accepted_steps)/(accepted_steps+rejected_steps));

}





int main()
{
  srand(time(NULL));   /* initialize random seed */


  double *xfield;
  xfield = my_vector(2*x_size*y_size);

  for(int i=0; i<2*x_size*y_size; i++) { xfield[i] = 0.0 ; }  /* initialize xfield to zero at all lattice sites for real and imag parts*/


  FILE *outsample ;
  char filename[300];


  
 
  /* printf("%f\n", sqrt(initial_mean_density)); */
  
 
  /* pick a randomly chosen site and updates field values mc_total_steps_site times for the lattice  */


  mc_updated_state(xfield, mc_total_steps_site); 



  for(int number_sample=0; number_sample<total_samples; number_sample++)
    {
      mc_updated_state(xfield, new_mc_updates); /* Again updating the field new_mc_upadtes times for a randomly chosen site */
    
      
      /* storing data to files from monte carlo generated initialization for each sample */
      sprintf(filename, "mc_field/mc_state_field_2d_%d", number_sample); 
      outsample = fopen(filename, "w");
      for(int i=0; i<2*x_size*y_size; i++)
	{
	  fprintf(outsample,"%le\n", xfield[i]);
	}
      fclose(outsample);
      
      

    }

  free(xfield);
  return 0;

}