myprojects-hg-reborn

Fork

// rk4_fixed.cpp 

/*
  Function to advance set of coupled first-order o.d.e.s by single step
  using fixed step-length fourth-order Runge-Kutta scheme

     x        ... independent variable
     y        ... array of dependent variables 
     h        ... fixed step-length

  Requires right-hand side routine
  
     void rhs_eval (double x, Array<double,1> y, Array<double,1>& dydx)
        
  which evaluates derivatives of y (w.r.t. x) in array dydx
*/




//////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////

/// This is a simple code implementing a RK fixed step 4th order method
/// to solve the harmonic oscillator equation (2nd order equation).


/////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////



#include <stdio.h>
#include <math.h>
#include <iostream>            
#include <cmath>
#include <vector>
#include <algorithm>
#include <complex>
#include <fstream>
#include "mylib.h"


#include <blitz/array.h>
/* Define complex double data type */


/* Define complex double data type */
typedef complex<double> dcomp;

using namespace blitz;
/*
void calc_cent_site(double Jt, double U, Array<complex<double>,1> inistate, vector<double> & Cinteraction,vector<double> & Ckin)
{int nsites=inistate.extent(0);
double onehalf=0.5;
int middle=(nsites-1)/2;  // remember that the first element in an array is arr(0) and not arr(1)

double number = -onehalf*U*abs(inistate(middle))*abs(inistate(middle))
*abs(inistate(middle))*abs(inistate(middle));
Cinteraction.push_back(number);

number=Jt*abs(inistate(middle)-inistate(middle-1))
*abs(inistate(middle)-inistate(middle-1))
+Jt*abs(inistate(middle+1)-inistate(middle))
*abs(inistate(middle+1)-inistate(middle))
;

Ckin.push_back(number);



}





void calc_lattice(double Jt, double U, Array<complex<double>,1> inistate,vector<double> & Tinteraction,vector<double> & Tkin)
{
int nsites=inistate.extent(0);
double onehalf=0.5;
double number1 =0;
double number2 =0 ;
for (int i=0;i<nsites;i++)
{

number1=number1-onehalf*U*abs(inistate(i))*abs(inistate(i))*abs(inistate(i))*abs(inistate(i));

// the previous expression accounts for the total interaction



if (i == 0)
{
number2=number2+Jt*abs(inistate(i)-inistate(nsites))*abs(inistate(i)-inistate(nsites));
}
else
{
number2=number2+Jt*abs(inistate(i)-inistate(i-1))*abs(inistate(i)-inistate(i-1));
}

  



}

Tinteraction.push_back(number1); 
Tkin.push_back(number2);

}


*/


int main ()
{
/*
double x=0;
 double h=0.01;
 


 Array<double,1> y(2); 
 Array<double,1> dydx(2);

 */
 
/*
// Parameters for the adaptive method

double acc=0.00002;
double S=10;
int maxrept=20;
double h_min=0.00000005;;
double h_max=0.0000025;;
int flag=1;
int rept;
double t_err;
*/

//Physical Parameters

 double x=0;
 double h=0.00000125;
 double U=-1.0;
double Jt=-1.0;
 double time=2.0;
        int nsites=51;
int nsave=200; //how many density profiles I want to save


double ratio=time/h/nsave;
cout<<"ratio is "<<ratio<<endl;
int every=rint(ratio);
cout<<"every is "<<every<<endl;
int niter=every*nsave;
cout<<"niter is "<<niter<<endl;
//Now a vector to use to save the population of the central well at different times
vector<double> savecen;
//vector<double> savetime;

int periodic=1; // 1 for periodic boundary conditions
                // and zero for non-periodic boundary conditions.
        
 Array<complex<double>,1> y(1); 
 Array<complex<double>,1> dydx(1);
 
// Now an array giving the initial state        
 Array<complex<double>,1> inistate(nsites),temp(nsites);
 
 cout<<"the number of sites is "<<nsites<<endl;
 Array<double,1> density(nsave*nsites);
 int middle=(nsites -1)/2;

 

 dcomp cn=dcomp(0,0);
 inistate=cn;
 cn=dcomp(sqrt(9.0),0);
/*
 for (int i=middle;i<(middle+1);i++)
 {
         inistate(i)=cn;
 }
 
 for (int i=0;i<nsites;i++)
 {
         cout<<"inistate is "<<inistate(i)<<endl;
 }
 
*/

// I will now read the initial state from a file and then normalize it to abs(cn)


vector<double> mydata;

string line;

  ifstream myfile ("./modeone.txt");

if (myfile.is_open())

{

while (!
myfile.eof() )

{   getline (myfile,line);

 cout << line << endl;


// ecco la riga! d è il tuo double
// double d = strtod( line.c_str() );

double d = strtod( line.c_str(), NULL );

 mydata.push_back(d);

}

myfile.close();


}


else cout <<
"unable to open file";

// now the file is stored in the mydata array
// and is normalized to one by default


 for (int i=0;i<nsites;i++)
{
inistate(i)=mydata[i]*cn;

}



double finpop[nsites];








double norm=0;
 for (int i=0;i<nsites;i++)
 {
         norm=norm+abs(inistate(i))*abs(inistate(i));
 finpop[i]=abs(inistate(i))*abs(inistate(i));


	 }
 
	 my_save("inipop.txt",nsites,  finpop);
	 cout << "the norm is "<<norm<<endl;






// Now I create an array saving the parameters of the simulation

vector<double> param;

param.push_back(Jt);
param.push_back(U);
param.push_back(h);
param.push_back(time);
 

                              
                            
 
int length=param.size();
my_save_vec("param.txt",length,  param);



 vector<double> Tinteraction; // vector containing the total interaction in the system 
                              // as a function of time

vector<double> Tkin; // vector containing the total kinetic in the system 
                              // as a function of time

vector<double> Cinteraction; // vector containing the  interaction in the central site system 
                              // as a function of time

vector<double> Ckin; // vector containing the kinetic interaction in the system 
                              // as a function of time

 int mycount =0;


// double onehalf =0.5;

while (x<=time)
 {
if (mycount % every == 0)
{
double cenpop=abs(inistate(middle))*abs(inistate(middle));
savecen.push_back(cenpop);


// Now I test my new functions

calc_cent_site( Jt,  U,  inistate,  Cinteraction, Ckin);

 calc_lattice( Jt, U, inistate, Tinteraction,Tkin) ;



/*
double number = -onehalf*U*abs(inistate(middle))*abs(inistate(middle))
*abs(inistate(middle))*abs(inistate(middle));
Cinteraction.push_back(number);

number=Jt*abs(inistate(middle)-inistate(middle-1))
*abs(inistate(middle)-inistate(middle-1))
+Jt*abs(inistate(middle+1)-inistate(middle))
*abs(inistate(middle+1)-inistate(middle))
;

Ckin.push_back(number);

*/

/*
double number1 =0;
double number2 =0 ;
for (int i=0;i<nsites;i++)
{

number1=number1-onehalf*U*abs(inistate(i))*abs(inistate(i))*abs(inistate(i))*abs(inistate(i));

 


if (i == 0)
{
number2=number2+Jt*abs(inistate(i)-inistate(nsites))*abs(inistate(i)-inistate(nsites));
}
else
{
number2=number2+Jt*abs(inistate(i)-inistate(i-1))*abs(inistate(i)-inistate(i-1));
}

  



}

Tinteraction.push_back(number1); 
Tkin.push_back(number2);
*/

}
        

 for (int j=0;j<nsites;j++)
{
        
         y(0)=inistate(j);




 rk4_fixed (x, y,  rhs_eval2,  h, j, inistate, U,  Jt, periodic);


// rk4_fixed (x, y,  rhs_eval,  h, j, inistate, U,  Jt, periodic);

temp(j)=y(0);

 }       
inistate=temp;   

// now I use inistate (containing the updated state of the system) 
// to work out the various interaction and energy parameters




// 
 mycount=mycount+1;
 }
 
 
 cout<<"time is "<<x<<endl;
 for (int i=0;i<nsites;i++)
 {
cout<<"the final state is "<<inistate(i)<<endl;
 }       
 
 
 norm=0;
 for (int i=0;i<nsites;i++)
 {
         norm=norm+abs(inistate(i))*abs(inistate(i));
 finpop[i]=abs(inistate(i))*abs(inistate(i));
	 }
 
	 my_save("finpop.txt",nsites,  finpop);
	 cout << "the norm is "<<norm<<endl;

length=savecen.size();
 
my_save_vec("cenpopn.txt",length,  savecen);

length=Tkin.size();
my_save_vec("Kin.txt",length,  Tkin);
my_save_vec("Utot.txt",length,  Tinteraction);

my_save_vec("CKin.txt",length,  Ckin);
my_save_vec("CUtot.txt",length,  Cinteraction);

return 0;
}