/*
 SimTest -- test out sim classes
 Barry Isralewitz   2002-Jan-18
 */
 

import java.awt.*;
import java.applet.*;
import java.util.*;
import java.text.*;
import langevinSimResult;
import cubbyHole;


public class langevinSim  {

    //public simResult result; 

    //only want one Random seed for all instances
    private static long seed;
   
    public int   timestep;
    private double param1;
    private double   Pi=3.1415926535897932e0;
    public int     C_units;
    private double   dt, dx; 
    public double D_r, D_p, R,R_rotor,R_arg210;
    public double lambda, eps_stator,R_pr;
    public double pH_A, pH_B, pK_A; 
    public double phi_A, phi_B, psi;	
    public double tau;      

    private double theta;
    private int pie;
    private int state; 
     //holds the information of two booleans
    private double theta_start;
    public double Rotor_step;
    public double theta_3;
    public double theta_1;
    public double theta_2;

    private double dist_L;
    private double El_pot_R0;
  
    public int step;
        //very hacky, set to N_steps to stop a  run...   	    
	//System.out.println("El_pot_R0=" + El_pot_R0);

    private NumberFormat myFormat = NumberFormat.getNumberInstance() ;

//constructor
  public langevinSim (int par1, double par2, int par3, int par4, long par5, int par6,double par7,double par8, double par9, double par10, double par11, double par12, double par13, double par14, double par15, double par16, double par17,double par18, double par19, double par20, double par21, double par22, double par23) {

	timestep = par1;
	theta    = par2;
	pie      = par3;
	state    = par4;
	seed     = par5;
	C_units  = par6; 
	dt       = par7;
	dx       = par8;
	D_r      = par9;
	D_p      = par10;
	R        = par11;
	R_rotor  = par12;
	R_arg210 = par13;
	lambda   = par14;
	eps_stator= par15;
	R_pr     = par16;
	pH_A     = par17;
	pH_B     = par18;
	pK_A     = par19;
	phi_A    = par20;
	phi_B    = par21;
	psi      = par22;
	tau      = par23;

	// Random my_random_numbers = new Random(seed);

    //// should be moved/copied to runSteps so can handle change in C_units
    Rotor_step=2.e0*Pi/((double)C_units);
	theta_3 =7.9e0/4.e0*Rotor_step/2.e0;
	theta_1=  Rotor_step/2.e0 - R/(2.e0*Pi*R_rotor)*2.e0*Pi;
	theta_2=  Rotor_step/2.e0 + R/(2.e0*Pi*R_rotor)*2.e0*Pi;

	//	System.out.println(Rotor_step+ "   "+ El_pot_R0+ " " +  R_arg210);
  }
    
    public langevinSimResult runSteps(int N_step) {

    //System.out.println("Starting sim, timestep= "+timestep+" theta="+ theta );
    //calc, in case C_Units has been changed by user
    //put back in when know how to deal with displacement too large error.
	Rotor_step=2.e0*Pi/((double)C_units);
	theta_3 =7.9e0/4.e0*Rotor_step/2.e0;
	theta_1=  Rotor_step/2.e0 - R/(2.e0*Pi*R_rotor)*2.e0*Pi;
	theta_2=  Rotor_step/2.e0 + R/(2.e0*Pi*R_rotor)*2.e0*Pi;
     
	theta_start=theta;
	double k_Ain_0  = 0.6e0*Math.pow(10.e0,-pH_A)*4.e0*R*D_p*Math.exp( phi_A );
	double k_Aout_0 = 0.6e0*Math.pow(10.e0,-pK_A)*4.e0*R*D_p*Math.exp((phi_A-phi_B-psi)/2.e0);
	double k_Bin_0  = 0.6e0*Math.pow(10.e0,-pH_B)*4.e0*R*D_p*Math.exp( phi_B );
	double k_Bout_0 = 0.6e0*Math.pow(10.e0,-pK_A)*4.e0*R*D_p*Math.exp((psi-phi_A+phi_B)/2.e0);
    //System.out.println("In runSteps, theta= "+theta+" pK_A ="+pK_A);

	//------Compute Electrostatic potential at the stator ends      
	double dist_L = Math.sqrt(R_arg210*R_arg210+R_pr*R_pr - 2.e0*R_arg210*R_pr*Math.cos(theta_1 ));
	double El_pot_R0= 56.23e0/dist_L/eps_stator*Math.exp(-lambda*dist_L);

    //starts a new Random each call (seed is incremented each call as well),
   // if costly, can make a static object for this...
	 Random my_random_numbers = new Random(seed);

	//-----zero sum_'s for statistics



	double sum_E=0.e0;
	double sum_R=0.e0;
	double sum_F=0.e0;
	double sum_L=0.e0;

	//-----now iterate 

	int sum = 0;
	double theta_old=0.e0;
	double asp61R= Rotor_step/2.e0;
	double asp61L= -Rotor_step/2.e0;
	double dist_R=0.e0;
	double psi_E_R=0.e0;
	double psi_E_L=0.e0;
	double k_Ain=0.e0;
	double k_Aout=0.e0;
	double k_Bin=0.e0;
	double k_Bout=0.e0;
	double K_ER = 0.e0;
	double K_LF = 0.e0;
	double K_RE = 0.e0;
	double K_FL = 0.e0;
	double K_EL = 0.e0;
	double K_RF = 0.e0;
	double K_LE = 0.e0;
	double K_FR = 0.e0;
	double[][] K = new double[4][4];
	int state_old=0;
	int asp61L_state = 0;
	int asp61R_state = 0;
	double angle_factor_R = 0.e0;
	double angle_factor_L = 0.e0;
	double F_eR = 0.e0;
	double F_eL = 0.e0;
	double F_e = 0.e0;
	int pie_old = 0;
	double prev=theta - Rotor_step*(double)pie;

	
      for (step=0; step < N_step; step++ ) {
            //System.out.println("Starting sim, timestep= "+timestep+" theta="+ theta );

	    theta_old = theta;

	    //--------check if displacement is larger than half of Rot_step (just in case)
	    if ((prev > (Rotor_step/2.e0)) | (prev<(-Rotor_step/2.e0))) { 
	    	String err = "displacement is larger than half of the rotor! displacement= " + myFormat.format(prev) ;
		 return new langevinSimResult(timestep,theta,pie,state,err);
	    };

	    //--------calculate theta_distance
	   
	    prev = theta - Rotor_step*(double)pie;
	    
	    //--------compute asp61L--arg210 and asp61R--arg210 distances
   
	    dist_R= Math.sqrt(R_arg210*R_arg210 + R_pr*R_pr - 2.e0*R_arg210*R_pr*Math.cos(Rotor_step/2.e0-prev));
	    dist_L= Math.sqrt(R_arg210*R_arg210 + R_pr*R_pr - 2.e0*R_arg210*R_pr*Math.cos(-Rotor_step/2.e0-prev));
	
	
	    //--------compute electrostatic potential:
	    
	    if ((Rotor_step/2.e0-prev)<theta_1) {
		psi_E_R=56.23e0/dist_R/eps_stator*Math.exp(-lambda*dist_R)-El_pot_R0;
	    } else {
		psi_E_R=0.e0; 
	    }
	    if ((Rotor_step/2.e0+prev)<theta_1) {
		psi_E_L=56.23e0/dist_L/eps_stator*Math.exp(-lambda*dist_L)-El_pot_R0;
	    } else {
		psi_E_L=0.e0;  
	    }
	    //--------compute transition rates (theta dependent)

	    k_Ain  = k_Ain_0;
	    k_Aout = k_Aout_0*Math.exp(psi_E_R);
	    k_Bin  = k_Bin_0;
	    k_Bout = k_Bout_0*Math.exp(psi_E_L);
	    
         
	    K_ER = k_Ain;
	    K_LF = k_Ain;
	    K_RE = k_Aout;
	    K_FL = k_Aout;
	    K_EL = k_Bin;
	    K_RF = k_Bin;
	    K_LE = k_Bout;
	    K_FR = k_Bout;
	    
	    K[0][0] = - (K_ER + K_EL);
	    K[0][1] = K_RE;
	    K[0][2] = 0.e0;
	    K[0][3] = K_LE;
	    K[1][0] = K_ER;
	    K[1][1] = - (K_RE + K_RF);
	    K[1][2] = K_FR;
	    K[1][3] = 0.e0;
	    K[2][0] = 0.e0;
	    K[2][1] = K_RF;
	    K[2][2] =  - (K_FR + K_FL);
	    K[2][3] = K_LF;
	    K[3][0] = K_EL;
	    K[3][1] = 0.e0;
	    K[3][2] = K_FL;
	    K[3][3] =  - (K_LE + K_LF);
	    //--------save old system state         
	    state_old = state;
	     //--------determine protonation states of each asp61 within the stator
	    if (state==1) {
		asp61L_state = 0;
		asp61R_state = 0;
	    } else if (state==2) {
		asp61L_state = 0;
		asp61R_state = 1;
	    } else if (state==3) {
		asp61L_state = 1;
		asp61R_state = 1;
	    } else if (state==4) {
		asp61L_state = 1;
		asp61R_state = 0;
	    }

	    //--------Deternime what will happen with asp61R:

            if (asp61R_state == 0) {
		if (k_Ain*dt > my_random_numbers.nextDouble() ) {
		    asp61R_state = 1; 
		}
	    } else {
		if (k_Aout*dt > my_random_numbers.nextDouble() ) {
		    asp61R_state = 0; 
		}
            }
	    //--------Deternime what will happen with asp61L:

            if (asp61L_state==0) {
		if (k_Bin*dt > my_random_numbers.nextDouble() ) {
		    asp61L_state = 1; 
		}
            } else {
		if (k_Bout*dt >  my_random_numbers.nextDouble() ) {
		    asp61L_state = 0; 
		}
            }    
	    //--------Determine state of the system:

	    if ((asp61L_state == 0) && (asp61R_state == 0)) {
		state = 1; }
	    if ((asp61L_state == 0) && (asp61R_state == 1)) {
		state = 2; }
	    if ((asp61L_state == 1) && (asp61R_state == 1)) {
		state = 3; }
	    if ((asp61L_state == 1) && (asp61R_state == 0)) {
		state = 4; }  

	    //--------write down the state of the system
	    // if (state_old != state) {
	    //	System.out.println("State has changed at step "+ step +"-> "+ state);
	    // }				     
	    //--------compute statistics 
	    if (state == 1) {
		sum_E=sum_E + 1.e0;
	    } else if (state == 2) {
		sum_R=sum_R + 1.e0;
	    } else if (state == 3) {
		sum_F=sum_F + 1.e0;
	    } else if (state == 4) {
		sum_L=sum_L + 1.e0;
	    } else {
		 String err = "state is not defined";
		 return new langevinSimResult(timestep,theta,pie,state,err);
	    }
	    //--------determine electrostatic forces:

	    angle_factor_R= R_arg210*R_pr*Math.sin(Rotor_step/2.e0-prev);
	    angle_factor_L=R_arg210*R_pr*Math.sin(-Rotor_step/2.e0-prev);

	    
	    if ((Rotor_step/2.e0-prev) < theta_1) {
		F_eR= 56.23/eps_stator/Math.pow(dist_R,3)*Math.exp(-lambda*dist_R)*(1.e0 + lambda*dist_R)*angle_factor_R;
	    } else {
		F_eR=0.e0;
	    }
	    if ((Rotor_step/2.e0+prev) < theta_1) {
		F_eL= 56.23/eps_stator/Math.pow(dist_L,3)*Math.exp(-lambda*dist_L)*(1.e0 + lambda*dist_L)*angle_factor_L;
	    } else {
		F_eL=0.e0;
	    }

	    if (state == 1) {
		F_e = F_eR + F_eL;
	    }
	    if (state == 2) {
		F_e = F_eL;
	    }
	    if (state == 3) { 
		F_e = 0.e0;
	    }
	    if (state == 4) {
		F_e = F_eR;
	    }
	    //--------Main iteration
	    //	    System.out.println(my_random_numbers.nextGaussian());
	    theta = theta_old + (F_e*D_r + tau*D_r + Math.sqrt(2.e0*D_r/dt)*my_random_numbers.nextGaussian()  )*dt; 
	    //--------Main output
	    
	    // System.out.println(timestep +" "+ theta +" "+ asp61L_state +" "+ asp61R_state);

	    //------DETERMINE pie_tracer ->
	    state_old = state;
	    pie_old = pie;

	    //--------If stator moves right...

	    if ((theta - Rotor_step*(double)pie) > Rotor_step/2.e0) { // stator moves right...
		//-----------change chemical state due to rotation
		if (state == 1) {
		    pie=pie_old;
		    theta = theta_old;
		    state=state_old;
		} else if (state == 2) {
		    pie=pie_old;
		    theta = theta_old;
		    state=state_old;
		} else if (state == 4) {
		    state = 2;        // Right
		    pie = pie + 1;
		} else if (state == 3) {
		    state = 3;        // Full
		    pie = pie + 1;
		} else {
		    String err="state of the system is not determined";
		    return new langevinSimResult(timestep,theta,pie,state,err);
		}            
		
		//--------If stator moves left...               
	    } else if ((theta - Rotor_step*(double)pie)<(-Rotor_step/2.e0)) {// stator moves left
		//-----------change chemical state due to rotation
		if (state == 1) {
		    pie=pie_old;
		    theta = theta_old;
		    state=state_old;
		} else if (state == 2) {
		    pie = pie - 1;
		    state = 4;        // Left
		} else if (state == 4) {
		    pie=pie_old;
		    theta = theta_old;
		    state=state_old;
		} else if (state == 3) {
		    pie = pie - 1;              
		    state = 3;        // Full
		} else {
		    String err = "state of the system is not determined";
		    return new langevinSimResult(timestep,theta,pie,state,err);
		}      
		//--------displacement is not sufficient to change the chemical state
	    } else {
		pie = pie;
	    }
	   
            timestep = timestep + 1;
	   
      }
   //increment seed so next run uses different value
	seed += 1;

		//System.out.println("Average rotation speed (Hz)="+ (theta -theta_start)/(dt*(double)N_step*2*Pi) +"; final theta="+ theta );
		//System.out.println("Leaving sim, timestep= "+timestep+" final theta="+ theta );
     return new langevinSimResult(timestep,theta,pie,state,"");
    }


//not needed, might use as a debugging check, but 

    public void printState() {
	System.out.println("timestep= " + timestep+ " theta= "+  theta +" pie=" + pie+ " state=" + state);
	
    } 
    
   public void doReset (int t, double th, int p, int s) {  
         timestep = t;
         theta = th;
         pie = p;
         state =s ;
   }
	
}