/* author: V. P. Singh */ #include #include #include #define i2d(x_site, y_site, x_length) (x_site + (y_site)*(x_length)) #define i2d_re(x_site, y_site, x_length) (2*x_site + 2*(y_site)*(x_length)) #define i2d_im(x_site, y_site, x_length) (2*x_site + (2*(y_site)*(x_length)) + 1) #define pn(x_site, x_length) (((x_length) + (x_site)%(x_length))%(x_length)) #define i2d_right_re(x_site, y_site, x_length, y_length) (2*pn((x_site+1), x_length) + 2*(y_site)*(x_length)) #define i2d_left_re(x_site, y_site, x_length, y_length) (2*pn((x_site-1), x_length) + 2*(y_site)*(x_length)) #define i2d_up_re(x_site, y_site, x_length, y_length) (2*x_site + 2*pn((y_site+1), y_length)*(x_length)) #define i2d_down_re(x_site, y_site, x_length, y_length) (2*x_site + 2*pn((y_site-1), y_length)*(x_length)) #define i2d_right_im(x_site, y_site, x_length, y_length) (2*pn((x_site+1), x_length) + (2*(y_site)*(x_length)) + 1) #define i2d_left_im(x_site, y_site, x_length, y_length) (2*pn((x_site-1), x_length) + (2*(y_site)*(x_length)) + 1) #define i2d_up_im(x_site, y_site, x_length, y_length) (2*x_site + (2*pn((y_site+1), y_length)*(x_length)) + 1) #define i2d_down_im(x_site, y_site, x_length, y_length) (2*x_site + (2*pn((y_site-1), y_length)*(x_length)) + 1) #define i2d_right_diag_re(x_site, y_site, x_length, y_length)(2*pn((x_site+1), x_length) + 2*pn((y_site+1), y_length)*(x_length)) #define i2d_right_diag_im(x_site, y_site, x_length, y_length)(2*pn((x_site+1), x_length) + (2*pn((y_site+1), y_length)*(x_length)) +1) #define i3d(x_site, y_site, z_site, x_length, y_length) (x_site + (y_site)*(x_length) + (z_site)*(x_length)*(y_length)) #define i2d_s2(x_site, y_site, x_length, y_length) (x_site + (y_site)*(x_length) + (x_length)*(y_length) ) #define i2d_s2_re(x_site, y_site, x_length, y_length) (2*x_site + 2*(y_site)*(x_length) + 2*(x_length)*(y_length) ) #define i2d_s2_im(x_site, y_site, x_length, y_length) (2*x_site + 2*(y_site)*(x_length) + 2*(x_length)*(y_length) + 1) #define i2d_s2_right_re(x_site, y_site, x_length, y_length) (2*pn((x_site+1), x_length) + 2*(y_site)*(x_length) + 2*(x_length)*(y_length)) #define i2d_s2_left_re(x_site, y_site, x_length, y_length) (2*pn((x_site-1), x_length) + 2*(y_site)*(x_length) + 2*(x_length)*(y_length)) #define i2d_s2_up_re(x_site, y_site, x_length, y_length) (2*x_site + 2*pn((y_site+1), y_length)*(x_length) + 2*(x_length)*(y_length) ) #define i2d_s2_down_re(x_site, y_site, x_length, y_length) (2*x_site + 2*pn((y_site-1), y_length)*(x_length) + 2*(x_length)*(y_length) ) #define i2d_s2_right_im(x_site, y_site, x_length, y_length) (2*pn((x_site+1), x_length) + (2*(y_site)*(x_length)) + 2*(x_length)*(y_length) + 1) #define i2d_s2_left_im(x_site, y_site, x_length, y_length) (2*pn((x_site-1), x_length) + (2*(y_site)*(x_length)) + 2*(x_length)*(y_length)+ 1) #define i2d_s2_up_im(x_site, y_site, x_length, y_length) (2*x_site + (2*pn((y_site+1), y_length)*(x_length)) + 2*(x_length)*(y_length)+ 1) #define i2d_s2_down_im(x_site, y_site, x_length, y_length) (2*x_site + (2*pn((y_site-1), y_length)*(x_length)) + 2*(x_length)*(y_length) + 1) #define i2d_s2_right_diag_re(x_site, y_site, x_length, y_length)(2*pn((x_site+1), x_length) + 2*pn((y_site+1), y_length)*(x_length) + 2*(x_length)*(y_length) ) #define i2d_s2_right_diag_im(x_site, y_site, x_length, y_length)(2*pn((x_site+1), x_length) + 2*pn((y_site+1), y_length)*(x_length) + 2*(x_length)*(y_length) +1 ) #define JMAX 17 #define JMAX1 30 #define JMAXP (JMAX+1) #define K 4 #define PI (double)3.141592653589793 #define EPS 1e-4 #define EPS1 1e-5 #define EPS2 1e-5 #define FUNC(x) ((*func)(x)) #define FUNC1(x) ((*func1)(x)) #define SAFETY 0.9 #define PGROW -0.2 #define PSHRNK -0.25 #define ERRCON 1.89e-4 #define MAXBIT 30 #define MAXDIM 6 static double xsav; static double (*nrfunc)(double,double); static char buf[12]; #define FREE_ARG char* char *itoa1(int i); char* random_state; double* my_vector(int array_size); void nrerror(const char error_text[]); double *vector(int n); void free_vector(double *v,int n); void rkqs(double y[], double dydx[], int n, double *x, double htry, double eps, double yscal[], double *hbid, double *hnext, void (*dervis)(double,double[],double[])); double* my_vector(int array_size) { return (double*)malloc(size_t (array_size*sizeof(double))); } double *vector(int n) { double *v; v=(double *)malloc((size_t) ((n+1)*sizeof(double))); return v+1; } void free_vector(double *v, int n) { free((FREE_ARG) (v-1)); } void nrerror(const char error_text[]) { fprintf(stderr,"Numerical Recipes run-time error...\n"); fprintf(stderr,"%s\n",error_text); fprintf(stderr,"...now exiting to system...\n"); exit(1); } void rkqs(double y[], double dydx[], int n, double *x, double htry, double eps, double yscal[], double *hbid, double *hnext, void (*dervis)(double,double[],double[])) { void rkck(double y[], double dydx[], int n, double x, double h, double yout[], double yerr[], void (*dervis)(double, double[], double[])); int i; double errmax,h, htemp, xnew, *yerr, *ytemp1; yerr=vector(n+1); ytemp1=vector(n+1); h=htry; while(true){ rkck(y,dydx,n,*x,h,ytemp1,yerr,dervis); errmax=0.; for(i=0;i<=n;i++) errmax=fmax(errmax,fabs(yerr[i]/yscal[i])); errmax/=eps; if(errmax<=1.0) break; htemp=SAFETY*h*pow(errmax,PSHRNK); h=(h>=0. ? fmax(htemp,0.1*h) : fmin(htemp,0.1*h)); xnew=(*x)+h; if (xnew==*x) nrerror("stepsize underflow in rkqs"); } if (errmax>ERRCON) *hnext=SAFETY*h*pow(errmax,PGROW); else *hnext=5.*h; *x+=(*hbid=h); for(i=0; i<=n; i++) y[i]=ytemp1[i]; free_vector(ytemp1,n+1); free_vector(yerr,n+1); } void rkck(double y[], double dydx[], int n, double x, double h, double yout[], double yerr[], void(*dervis)(double, double[], double [])) { int i; static double a2=0.2, a3=0.3,a4=0.6,a5=1.0,a6=0.875,b21=0.2,b31=3.0/40.0, b32=9./40., b41=0.3, b42=-0.9, b43=1.2, b51=-11./54., b52=2.5, b53=-70./27. , b54=35./27. , b61=1631./55296. , b62=175./512. , b63=575./13824., b64=44275./110592., b65=253./4096., c1=37./378.,c3=250./621., c4=125./594., c6=512./1771., dc5=-277./14336.; double dc1=c1-2825./27648., dc3=c3-18575./48384., dc4=c4-13525./55296., dc6=c6-0.25; double *ak2, *ak3, *ak4, *ak5, *ak6, *ytemp; ak2=vector(n+1); ak3=vector(n+1); ak4=vector(n+1); ak5=vector(n+1); ak6=vector(n+1); ytemp=vector(n+1); for (i=0; i<=n; i++) ytemp[i]=y[i]+b21*h*dydx[i]; (*dervis)(x+a2*h,ytemp,ak2); for (i=0; i<=n; i++) ytemp[i]=y[i]+h*(b31*dydx[i]+b32*ak2[i]); (*dervis)(x+a3*h,ytemp,ak3); for (i=0; i<=n; i++) ytemp[i]=y[i]+h*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]); (*dervis)(x+a4*h,ytemp,ak4); for (i=0; i<=n; i++) ytemp[i]=y[i]+h*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]); (*dervis)(x+a5*h,ytemp,ak5); for (i=0; i<=n; i++) ytemp[i]=y[i]+h*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]); (*dervis)(x+a6*h,ytemp,ak6); for (i=0; i<=n; i++) yout[i]=y[i]+h*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]); for (i=0; i<=n; i++) yerr[i]=h*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]); free_vector(ytemp,n); free_vector(ak6,n); free_vector(ak5,n); free_vector(ak4,n); free_vector(ak3,n); free_vector(ak2,n); } double gasdev(void){ static int iset=0; static double gset; double fac,rsq,v1,v2; if (iset == 0) { do { v1=2.0*random()/(RAND_MAX+1.0)-1.0; v2=2.0*random()/(RAND_MAX+1.0)-1.0; rsq=v1*v1+v2*v2; } while (rsq >= 1.0 || rsq == 0.0); fac=sqrt(-2.0*log(rsq)/rsq); gset=v1*fac; iset=1; return v2*fac; } else { iset=0; return gset; } } double trap_2d(int x_site, int y_site, int x_size, int y_size, double freq_x_Hz, double freq_y_Hz, double ener_factor) { return ener_factor*( pow( freq_x_Hz*(x_site - x_size/2.), 2) + pow( freq_y_Hz*(y_site - y_size/2.), 2) ) ; }