Skip to content

Instantly share code, notes, and snippets.

@ktchernov

ktchernov/gist:d553549b74994f5fd129

Created November 27, 2014 22:49
Show Gist options
  • Save ktchernov/d553549b74994f5fd129 to your computer and use it in GitHub Desktop.
Save ktchernov/d553549b74994f5fd129 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
gistfile1.c
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
#include<time.h>
#include "mpi.h"
#define PI 3.1415926535897931
#define OBSTYP 0
#define OBSM 1
#define dist2(x,x1,y,y1,z,z1) ((((x)-(x1))*((x)-(x1)))+(((y)-(y1))*((y)-(y1)))+(((z)-(z1))*((z)-(z1))))
#define crossx(x,y,z,x1,y1,z1) (((y)*(z1))-((z)*(y1)))
#define crossy(x,y,z,x1,y1,z1) (((z)*(x1))-((x)*(z1)))
#define crossz(x,y,z,x1,y1,z1) (((x)*(y1))-((y)*(x1)))
#define u11(ux,uy,uz,c,s) ((c)+(ux)*(ux)*((1.0)-(c)))
#define u12(ux,uy,uz,c,s) (((ux)*(uy)*((1.0)-(c)))-((uz)*(s)))
#define u13(ux,uy,uz,c,s) (((ux)*(uz)*((1.0)-(c)))+((uy)*(s)))
#define u21(ux,uy,uz,c,s) (((ux)*(uy)*((1.0)-(c)))+((uz)*(s)))
#define u22(ux,uy,uz,c,s) ((c)+(uy)*(uy)*((1.0)-(c)))
#define u23(ux,uy,uz,c,s) (((uy)*(uz)*((1.0)-(c)))-((ux)*(s)))
#define u31(ux,uy,uz,c,s) (((ux)*(uz)*((1.0)-(c)))-((uy)*(s)))
#define u32(ux,uy,uz,c,s) (((uy)*(uz)*((1.0)-(c)))+((ux)*(s)))
#define u33(ux,uy,uz,c,s) ((c)+(uz)*(uz)*((1.0)-(c)))
struct aftr
{
unsigned long int n,plane,pn,n1,n2,typ;
double xreal,yreal,zreal,x0,y0,z0,r;
double vx,vy,vz,ax,ay,az,fx,fy,fz,PE,hl,bn;
double round,body,cenben,kbd,fl;
double nx,ny,nz,nx0,ny0,nz0,tx,ty,tz,cangle1;
double cx,cy,cz,epx,emx,epy,emy,epz,emz;
double xa,ya,za,mass,rad,rad2,D;
struct aftr *flagella[2],*helix[2];
};
struct CM
{
double xreal,yreal,zreal,xc0,yc0,zc0;
double vx,vy,vz,fx,fy,fz,cx,cy,cz;
};
struct cnline
{
double xreal,yreal,zreal;
double dirx,diry,dirz;
double radius2,radius;
double cx,cy,cz;
};
struct sol
{
unsigned long int n,in,pn;
double x,y,z,vx,vy,vz;
double xa,ya,za;
};
struct wallp
{
unsigned long int n,in,pn;
double x,y,z,vx,vy,vz,vrx,vry,vrz;
double xa,ya,za;
};
struct redcell
{
unsigned long int id,ip,partnr;
double xcm,ycm,zcm,rx,ry,rz;
double ex,ey,ez;
double Vxcm,Vycm,Vzcm;
double wx,wy,wz;
double Momx,Momy,Momz,Amomx,Amomy,Amomz,coltim;
};
struct cellVCM
{
double Vx,Vy,Vz,Vnx,Vny,Vnz;
double S;
double xcm,ycm,zcm,crx,cry,crz;
double r11,r12,r13,r21,r22,r23,r31,r32,r33;
double xx,yy,zz;
};
/*****************************/
unsigned long int L,Lr1,length,seed,npoints,nflagella;
double Lr0;
double radius,k0,kbend,dtaf,TEMP,gamma1,wave,frequency,wavelength,ylc,zlc;
double Rrbc,Mrbc,Irbc;
/********************************************/
//////////////////////////////////////////
double starttimempi,endtimempi,totaltimempi;
int rank,numprocs,left,right,false,true,xer,xsl;
MPI_Comm comm1d;
MPI_Status status;
////////////////////////////////////////////////////
char foutpoint[256] = "davod";
char foutpoint1[256] = "davoda";
char fname[256] = "davoda";
/********************************************/
unsigned long int Naf,tsim,NAF,classwidth,helix,flalen,Nsol,Nsolp,Nsolw,Nrbc,Nrbct,Nrbcp;
unsigned long int *N_monocell,cellsize,nprcell,svent;
double angle,SRTEMP,tottime,sangle,eps,totdt,inv_Naf;
double R2CM,R2,KE,PE,dt;
double centerlength,totforx,totfory,totforz,nfx,nfy,nfz,totvx,totvy,totvz,tottorx,tottory,tottorz;
double vsol2,wv,fr,TNF1,TNF2,TPF,TNV1,TNV2,TPV,avgspcur;
double kbody,kround,kcenben,knew,kcomp,tlength,inv_tlength,kb1,fr1,kmi;
double tbf,taf,vaddx,vaddy,vaddz,totangx,totangy,totangz,inv_npoints,reducedbend;
/*********************************************/
struct aftr *at;
struct CM *cm;
struct cnline *cnt;
struct sol *ss;
struct wallp *ww;
struct redcell *rbc;
struct cellVCM *N_VCM;
/*********************************************/
FILE *fp;
time_t date;
char name[50];
/*********************************/
void init(void);
void afststraight(void);
void aftrvel1(void);
void lenat0(void);
void normal1(void);
void restart(void);
void restartb(void);
void transverse(void);
void rbcpos(void);
void solventposvel1(void);
void wallpos(void);
void streaming(double);
void snapshot(void);
void snapshot1(void);
void radiuscal(void);
void bendbodyforcet0(void);
void insideaftr8(unsigned long int);
void bounceback(int,unsigned long int,double);
void bouncebackrbc(int,unsigned long int,int,double);
void ND(void);
void MDRBC(void);
double gaussianrandom(void);
void calculateforces(void);
void calculateroundspring(void);
void bodyspring(void);
void cenbenspring(void);
void bendbodyforce2(void);
void forcebend11(void);
void collisionAT(void);
void average(void);
void results(void);
void flalink(struct aftr*,struct aftr*);
void hellink(struct aftr*,struct aftr*);
void velupdate(double,double,double,double,double,double,unsigned long int,unsigned long int);
void velupdate2(unsigned long int,unsigned long int);
double rayleighrandom(void);
/*****************************************/
int main(int argc, char *argv[])
{
unsigned long int i,n;
double clt,mass;
//////////////////////parallel//////////////////////
// double starttimempi,endtimempi,totaltimempi;
// int rank,numprocs,left,right,false,true;
false=0;
true=1;
////////////////////////////////////////////////////
//////////////////////parallel/////////////////////
/// MPI_Comm comm1d;
// comm1d=MPI_COMM_WORLD;
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
MPI_Cart_create(MPI_COMM_WORLD,1,&numprocs,&true,true,&comm1d);
MPI_Comm_rank(comm1d,&rank);
MPI_Cart_shift(comm1d,0,1,&left,&right);
///////////////////////////////////////////////////
// printf("\n rank=%i numprocs=%i \n",rank, numprocs);
// starttimempi=MPI_Wtime();
/////////////////////////////////////////////////
/* if(argc!=20)
{
printf("PB arguments != 19\n");
exit(1);
}
*/
Lr1=18;
L=120;
seed=0;
length=26;
radius=1.0;
npoints=10;
nflagella=0;
k0=3000000.0;
kbend=50.0;
dtaf=0.0001;
TEMP=1.0;
tottime=100000.0;
gamma1=1.0;
wave=1;
helix=1;
dt=0.01;
nprcell=10;
svent=1;
frequency=-0.005;
wavelength=2.25;
if(wave>2 || wave<0)
{
printf("african trypanosome change wave\n");
exit(1);
}
if(helix>1 || helix<0)
{
printf("what change helix\n");
exit(1);
}
Lr0=(Lr1-2.5)/2.0;
zlc=(double)Lr1/2.0;
ylc=(double)Lr1/2.0;
xsl=rank*L/numprocs;
xer=(rank+1)*L/numprocs;
if(L%numprocs != 0 )
{
exit(1);
}
cellsize=1;
centerlength=1.;
NAF=1;
Nrbc=10000;
srand48(seed);
SRTEMP=sqrt(TEMP);
Naf=(length*(npoints));
flalen=length;
Nsol=nprcell*L*Lr1*Lr1;
Nsolp=nprcell*Lr1*Lr1*L/numprocs;
mass=1.;
inv_Naf=1./(double)Naf;
inv_npoints=1./(double)npoints;
angle=2.*PI/(double)npoints;
wv=wavelength*2.*PI/(centerlength*((double)length-1.));
fr=-2.*PI*frequency;
reducedbend=0.8;
knew=10.;
kb1=180.;
kcenben=1.;
kbody=1.;
kround=1.;
kcomp=0;
kmi=180.;
avgspcur=0;
if(rank==0){
fp=fopen("info.txt","w");
date=time(NULL);
fprintf(fp,"Created by %s at %s",argv[0],ctime(&date));
fprintf(fp,"knew=%le\tkb1=%le\tkcenben=%le\tkbody=%le\tkround=%le\tkcomp=%le\tkmi=%le\n",knew,kb1,kcenben,kbody,kround,kcomp,kmi);
fclose(fp);
}
sangle=sin(angle*0.5);
eps=10.;
totdt=classwidth=0;
init();
afststraight();
aftrvel1();
lenat0();
bendbodyforcet0();
tsim=0;
if(tsim != 0)
{
restart();
Rrbc=1.0;
Mrbc=100*nprcell*4*PI*Rrbc*Rrbc*Rrbc/3;
Irbc=2*Mrbc*Rrbc*Rrbc/5;
// restartb();
clt=totdt+dt;
}
// tsim++;
// snapshot1();
// tsim--;
normal1();
//transverse();
if(tsim == 0)
{
rbcpos();
solventposvel1();
streaming(dt);
snapshot();
//radiuscal();
//bendbodyforcet0();
tsim=totdt=classwidth=0;
clt=dt;
//tlength=sqrt(dist2((at+0)->xreal,(at+length-1)->xreal,(at+0)->yreal,(at+length-1)->yreal,(at+0)->zreal,(at+length-1)->zreal));
//inv_tlength=1./tlength;
}
R2CM=R2=KE=PE=totforx=totfory=totforz=totvx=totvy=totvz=TNF1=TNF2=TPF=TNV1=TNV2=TPV=tottorx=tottory=tottorz=taf=0;
totangx=totangy=totangz=0;
starttimempi=MPI_Wtime();
//exit(1);
do
{
tsim++;
ND();
// dt=dtaf;
// MDRBC();
if(totdt>clt)
{
dt=0.01;
//printf("rankkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk=%d tsim = %i\n",rank,tsim);
MDRBC();
//printf("rankkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk22222222222222=%d tsim = %i\n",rank,tsim);
normal1();
///transverse();
vaddx=vaddy=vaddz=0;
streaming(dt);
if(OBSM==1) MDRBC();
if(numprocs > 1){
MPI_Allreduce(&vaddx,&vaddx,1,MPI_DOUBLE,MPI_SUM,comm1d);
MPI_Allreduce(&vaddy,&vaddy,1,MPI_DOUBLE,MPI_SUM,comm1d);
MPI_Allreduce(&vaddz,&vaddz,1,MPI_DOUBLE,MPI_SUM,comm1d);
}
vaddx*=inv_Naf;
vaddy*=inv_Naf;
vaddz*=inv_Naf;
collisionAT();
clt+=dt;
}
/// results();
totdt+=dtaf;
classwidth++;
if(tsim%1000==0)
{
/// average();
classwidth=0;
}
if(rank==0 && tsim%1000==0) printf("rank=%d tsim = %i\n",rank,tsim);
if(tsim%150000==0)
snapshot();
if(tsim%5000000==0)
snapshot1();
}
while(totdt<=tottime);
//////////////////////////////////////
endtimempi=MPI_Wtime();
totaltimempi=endtimempi-starttimempi;
printf("rank=%d totaltimempi = %e\n",rank,totaltimempi);
//////////////////////////////////////
MPI_Finalize();
/////////////////////////////////////////////
return 0;
}
void init()
{
unsigned long int mem,AT,SRD;
double ppb,c,alpha;
at=(struct aftr *)calloc(Naf,sizeof(struct aftr));
cm=(struct CM *)calloc(NAF,sizeof(struct CM));
cnt=(struct cnline *)calloc(flalen,sizeof(struct cnline));
ss=(struct sol *)calloc(Nsolp*2,sizeof(struct sol));
ww=(struct wallp *)calloc(Nsolp,sizeof(struct wallp));
rbc=(struct redcell *)calloc(Nrbc,sizeof(struct redcell));
N_VCM=(struct cellVCM *)calloc(Lr1*Lr1*L/numprocs,sizeof(struct cellVCM));
N_monocell=(unsigned long int *)calloc(Lr1*Lr1*L/numprocs,sizeof(unsigned long int));
mem=0;
mem+=Naf*sizeof(struct aftr)+NAF*sizeof(struct CM)+flalen*sizeof(struct cnline)+Nsolp*sizeof(struct sol);
mem+=L*L*L*(sizeof(struct cellVCM)+2*sizeof(unsigned long int))+3*sizeof(double);
ppb=(double)(Naf+Nsol)/(double)(L*L*L);
c=ppb/(ppb-1.+exp(-ppb));
AT=1;
if(rank==0){
fp=fopen("info.txt","a");
fprintf(fp,"L=%ld\nseed=%ld\nlenght=%ld\tradius=%le\tnpoints=%ld\tNaf=%ld\nk0=%le\tkbend=%le\t",L,seed,length,radius,npoints,Naf,k0,kbend);
fprintf(fp,"gamma1=%le\n",gamma1);
fprintf(fp,"dt=%le\tdtaf=%le\n",dt,dtaf);
if(SRD==1)
{
fprintf(fp,"Nsol=%ld\tppb=%le\nkineticviscosity=%le\tcollisionviscosity=%le\t",Nsol,ppb,((5.*c/(2.-cos(alpha)-cos(2.*alpha)))-1.)*TEMP*0.5*dt,((1-cos(alpha))/(18.*c))*(double)(cellsize*cellsize)/dt,(TEMP*dt*(ppb/(ppb-1.+exp(-ppb))-0.5))+((ppb-1.+exp(-ppb))/(ppb*12.*dt)));
fprintf(fp,"totalvicosity=%le\n",(((5.*c/(2.-cos(alpha)-cos(2.*alpha)))-1.)*TEMP*0.5*dt)+(((1-cos(alpha))/(18.*c))*(double)(cellsize*cellsize)/dt));
}
if(AT==1)
{
fprintf(fp,"Nsol=%ld\tppb=%le\nkineticviscosity=%le\tcollisionviscosity=%le\t",Nsol,ppb,((ppb/(ppb-(5./4.)))-0.5)*TEMP*dt,cellsize*cellsize/(24.*dt)*((ppb-(7./5.))/ppb));
fprintf(fp,"totalvicosity=%le\n",(((ppb/(ppb-(5./4.)))-0.5)*TEMP*dt)+(cellsize*cellsize/(24.*dt)*((ppb-(7./5.))/ppb)));
}
fprintf(fp,"wavelenght=%le\tfrequency=%le\tavgspontaneouscurvature=%le\n",2*PI/wv,fr,avgspcur);
fprintf(fp,"\n\n");
if(helix==0)
{
if(wave==0)
fprintf(fp,"Straight flagella with longitudnal wave passing along the flagella\n");
if(wave==1)
fprintf(fp,"Straight flagella with transverse wave passing along the flagella\n");
if(wave==2)
fprintf(fp,"Straight flagella with helical wave passing along the flagella\n");
}
if(helix==1)
{
if(wave==0)
fprintf(fp,"Helical flagella with longitudnal wave passing along the flagella\n");
if(wave==1)
fprintf(fp,"Helical flagella with transverse wave passing along the flagella\n");
if(wave==2)
fprintf(fp,"Helical flagella with helical wave passing along the flagella\n");
}
fprintf(fp,"\n\n");
fprintf(fp,"Memory allocation (MB) : %lf\n",(double)mem/(double)(1024*1024));
fclose(fp);
fp=fopen("angle.dat","w");
fclose(fp);
fp=fopen("diffusion.dat","w");
fclose(fp);
fp=fopen("for-vel.dat","w");
fclose(fp);
fp=fopen("fvpp.dat","w");
fclose(fp);
fp=fopen("tor.dat","w");
fclose(fp);
fp=fopen("sol-vel.dat","w");
fclose(fp);
}
}
void afststraight()
{
unsigned long int i,j,k,j1,j2,n,n1;
double x,y,z,r,phi,c,s,d,l0,p,l,t1,t2,t3,rad,r2,ccc;
double e,r1,ll,coangle,la;
struct aftr *a,*a1,*a2;
k=rad=0;
phi=l=0;
t2=t3=(double)Lr1*0.5;
t1=(double)L*0.065;
for(j=0;j<npoints;j++)
{
s=sin(phi);
c=cos(phi);
i=l=ccc=0;
do
{
a=(at+k);
a->n=k;
// r=(0.16*(double)i+0.15)*pow((0.3*(double)i+2.5),0.13)*exp(-0.0115*((double)i+0.8)*((double)i+0.8));
// r=1.75*(0.165*(double)i/1.5+0.15)*pow((0.3*(double)i/1.9+2.5),0.13)*exp(-0.0115*((double)i/1.5+0.8)*((double)i/1.5+0.8));
if(i==0){
r=0.5*radius*pow(((double)i+0.15),0.25)*exp(-0.006*((double)i+0.15)*((double)i+0.15));
}
else if(i<21){
// r=radius*pow(((double)i+0.15),0.25)*exp(-0.006*((double)i+0.15)*((double)i+0.15));
r=radius*pow(((double)i-1.+0.15),0.25)*exp(-0.006*((double)i-1.+0.15)*((double)i-1.+0.15));
}
else{
r=r*0.95;
}
a->r=r;
a->mass=1.;
a->kbd=kbend;
if(r>rad)
rad=r;
if(i==0)
{
l=0;
}
else
{
// r1=(0.16*(double)i+0.15)*pow((0.3*(double)i+2.5),0.13)*exp(-0.0115*((double)i+0.8)*((double)i+0.8));
// r1=1.75*(0.165*(double)(i-1)/1.5+0.15)*pow((0.3*(double)(i-1)/1.9+2.5),0.13)*exp(-0.0115*((double)(i-1)/1.5+0.8)*((double)(i-1)/1.5+0.8));
if(i==1){
//r=0.5*radius*pow(((double)i+0.15),0.25)*exp(-0.006*((double)i+0.15)*((double)i+0.15));
r1=0.5*radius*pow(((double)i-1.+0.15),0.25)*exp(-0.006*((double)i-1.+0.15)*((double)i-1.+0.15));
}
else if(i<22){
r1=radius*pow(((double)i-2.+0.15),0.25)*exp(-0.006*((double)i-2.+0.15)*((double)i-2.+0.15));
}
else{
r1=r1*0.95;
}
l+=sqrt(centerlength*centerlength-(r-r1)*(r-r1));
if(i==1)
ll=l;
if(i>1)
(at+k)->D=l-ll;
}
a->bn=r;
x=(at+k)->xreal=(at+k)->x0=t1+l;
y=(at+k)->yreal=(at+k)->y0=t2+r*s;
z=(at+k)->zreal=(at+k)->z0=t3+r*c;
if(i>0)
{
coangle=acos(1.*((at+k-1)->xreal-(at+k)->xreal))*180./PI;
ccc+=sqrt(dist2((at+k-1)->xreal,(at+k)->xreal,(at+k-1)->yreal,(at+k)->yreal,(at+k-1)->zreal,(at+k)->zreal));
// printf("%le\t%le\t%le\t%le\n",l,sqrt(dist2((at+k-1)->xreal,(at+k)->xreal,(at+k-1)->yreal,(at+k)->yreal,(at+k-1)->zreal,(at+k)->zreal)),ccc,180.-coangle);
}
/////////////////////////////////////
if(x >= xsl && x < xer){
a->typ=0;
}
else if ( x < xsl && (xsl-x) < (2*(xer-xsl))){
a->typ=1;
}
else if ( x >= xer && (x-xer) < (2*(xer-xsl))){
a->typ=1;
}
else{
a->typ=2;
}
// printf(" %i\t%i\t%le\t%i\t%i\n",rank,a->typ,a->xreal,xsl,xer);
///////////////////////////////////////////
if(i>0)
{
a1=(at+(k-1));
flalink(a1,a);
flalink(a,a1);
}
i++;
k++;
}
while(i<length);
phi+=angle;
}
la=0.0;
j=1;
a=at+2;
for(i=2;i<length;i++)
{
if(i<8)
{
a1=at+i+1;
a->hl=sqrt(dist2(a->xreal,a1->xreal,a->yreal,a1->yreal,a->zreal,a1->zreal));
hellink(a1,a);
hellink(a,a1);
if(a->n+182>260)
{
a->n1=a->n+182-260;
a->n2=a->n+78;
}
else
{
a->n1=a->n+182;
a->n2=a->n+78;
}
// a->D=sqrt(dist2((at+1)->xreal,a->xreal,(at+1)->yreal,a->yreal,(at+1)->zreal,a->zreal));
// la+=sqrt(dist2(a2->xreal,a->xreal,a2->yreal,a->yreal,a2->zreal,a->zreal));
a->D=la;
la+=sqrt(dist2(a->xreal,a1->xreal,a->yreal,a1->yreal,a->zreal,a1->zreal));
// printf("%ld\t%ld\t%ld\t%ld\t%le\n",a->n2,a->n,a->n1,a1->n,a->D);
}
if(i==9 || i==10 || i==11 || i==12 )
{
a1=(at+i+j*length);
//a->hl=sqrt(dist2(a->xreal,a1->xreal,a->yreal,a1->yreal,a->zreal,a1->zreal));
hellink(a1,a);
hellink(a,a1);
if(a->n+182>260)
{
a->n1=a->n+182-260;
a->n2=a->n+78;
}
else
{
a->n1=a->n+182;
a->n2=a->n+78;
}
a->D=la;
la+=sqrt(dist2(a->xreal,a1->xreal,a->yreal,a1->yreal,a->zreal,a1->zreal));
// a->D=sqrt(dist2((at+1)->xreal,a->xreal,(at+1)->yreal,a->yreal,(at+1)->zreal,a->zreal));
j++;
// printf("%ld\t%ld\t%ld\t%ld\t%le\n",a->n2,a->n,a->n1,a1->n,a->D);
}
if(i>12)
{
a1=(at+i+j*length);
//a->hl=sqrt(dist2(a->xreal,a1->xreal,a->yreal,a1->yreal,a->zreal,a1->zreal));
hellink(a1,a);
hellink(a,a1);
if(a->n+182>260)
{
a->n1=a->n+182-260;
a->n2=a->n+78;
}
else
{
a->n1=a->n+182;
a->n2=a->n+78;
}
// a->D=sqrt(dist2((at+4)->xreal,a->xreal,(at+4)->yreal,a->yreal,(at+4)->zreal,a->zreal));
a->D=la;
la+=sqrt(dist2(a->xreal,a1->xreal,a->yreal,a1->yreal,a->zreal,a1->zreal));
// printf("%ld\t%ld\t%ld\t%ld\t%le\n",a->n2,a->n,a->n1,a1->n,a->D);
}
if(i<length-1)
a=a1;
}
if(k!=Naf)
{
printf("pb %ld \t %ld\n",k,Naf);
exit(1);
}
radius=rad;
}
void aftrvel1()
{
unsigned long int i,j;
double Vx,Vy,Vz,xcm,ycm,zcm,v2,inv_temp,temp;
struct aftr *a;
Vx=Vy=Vz=xcm=ycm=zcm=0;
for(i=0;i<Naf;i++)
{
a=at+i;
a->vx=gaussianrandom();
a->vy=gaussianrandom();
a->vz=gaussianrandom();
Vx+=a->vx;
Vy+=a->vy;
Vz+=a->vz;
}
Vx*=1./(double)Naf;
Vy*=1./(double)Naf;
Vz*=1./(double)Naf;
v2=0;
for(i=0;i<Naf;i++)
{
a=at+i;
a->vx-=Vx;
a->vy-=Vy;
a->vz-=Vz;
v2+=a->vx*a->vx+a->vy*a->vy+a->vz*a->vz;
xcm+=a->xreal;
ycm+=a->yreal;
zcm+=a->zreal;
}
(cm+0)->xreal=(cm+0)->xc0=xcm*inv_Naf;
(cm+0)->yreal=(cm+0)->yc0=ycm*inv_Naf;
(cm+0)->zreal=(cm+0)->zc0=zcm*inv_Naf;
temp=v2/3.*inv_Naf;
inv_temp=sqrt(TEMP/temp);
v2=0;
for(i=0;i<Naf;i++)
{
a=at+i;
a->vx*=inv_temp;
a->vy*=inv_temp;
a->vz*=inv_temp;
v2+=a->vx*a->vx+a->vy*a->vy+a->vz*a->vz;
}
printf("TEMP=%le\n",1./3.*v2/(double)Naf);
}
void lenat0()
{
unsigned long int i,j,k;
double dam1x,dam1y,dam1z,fx,fy,fz,l0,d2,d,inv_d;
struct aftr *a1,*a2,*a3,*a4;
/******************************************************/
for(i=0;i<length;i++)
{
k=i;
j=0;
do
{
if(j==0)
{
a1=(at+k);
a2=(at+k+length*(j+1));
}
else
{
a1=a2;
if(j<npoints-1)
a2=(at+k+length*(j+1));
else
a2=(at+k);
}
dam1x=(a1->xreal-a2->xreal);
dam1y=(a1->yreal-a2->yreal);
dam1z=(a1->zreal-a2->zreal);
a1->round=sqrt(dam1x*dam1x+dam1y*dam1y+dam1z*dam1z);
j++;
}
while(j<npoints);
}
/******************************************************/
i=0;
do
{
for(j=0;j<length-1;j++)
{
a1=at+(i)*length+j;
a2=at+(i)*length+j+1;
dam1x=(a1->xreal-a2->xreal);
dam1y=(a1->yreal-a2->yreal);
dam1z=(a1->zreal-a2->zreal);
a1->body=sqrt(dist2(a1->xreal,a2->xreal,a1->yreal,a2->yreal,a1->zreal,a2->zreal));
}
i++;
}
while(i<npoints);
/*if(nflagella==1)
{
a1=at+length-1;
a2=at+Naf-1;
a1->body=sqrt(dist2(a1->xreal,a2->xreal,a1->yreal,a2->yreal,a1->zreal,a2->zreal));
}*/
/******************************************************/
for(i=0;i<length;i++)
{
k=i;
j=0;
if(i==0)
l0=2.*radius*0.5;
else
l0=2./(double)i*radius;
do
{
a1=(at+(k+j*length));
a2=(at+k+length*(j+5));
a1->cenben=sqrt(dist2(a1->xreal,a2->xreal,a1->yreal,a2->yreal,a1->zreal,a2->zreal));
a2->cenben=sqrt(dist2(a1->xreal,a2->xreal,a1->yreal,a2->yreal,a1->zreal,a2->zreal));
j++;
}
while(j<npoints-5);
}
}
void normal1()
{
unsigned long int i,j,k,m,ki=0,kj=0;
double d2,d,inv_d;
double dax,day,daz,x1;
struct aftr *a1;
struct cnline *c;
////////////////////////////////
static double *hud; //////////////
static double *hvd; //////////////
static double *had; //////////////
static double *hbd; //////////////
static int n3st = 0;
if(n3st == 0){
hud = (double *) malloc(sizeof(double)*(5000));///////////////
if(hud == NULL){
printf("\n err_message(2)");/////////////////////
exit(1);
}
hvd = (double *) malloc(sizeof(double)*(5000));/////////////
if(hvd == NULL) exit(1);////////////////
had = (double *) malloc(sizeof(double)*(5000));/////////////////
if(had == NULL) exit(1);/////////////////////
hbd = (double *) malloc(sizeof(double)*(5000));////////////////
if(hbd == NULL) exit(1);///////////////////
n3st++;
}
hud[ki++] = 0.0;
hvd[kj++] = 0.0;
had[0] = 0.0;
hbd[0] = 0.0;
/////////////////////////////////////////
for(i=0;i<Naf;i++)
{
a1=at+i;
if(a1->typ != 0)continue;
x1=a1->xreal;
if( x1 >= xer || x1 < xsl)
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRR1 \n");
printf(" rank=%i i=%i a1->typ=%i a1->xreal=%le xsl=%i xer=%i \n",rank,i,a1->typ,a1->xreal,xsl,xer);
exit(1);
}
hvd[kj++]=(double)i;
hvd[kj++]=a1->vx;
hvd[kj++]=a1->vy;
hvd[kj++]=a1->vz;
// printf("rank=%i i=%i ax=%le ay=%le az=%le vx=%le \n",rank,i,a1->ax,a1->ay,a1->az,a1->vx);
hud[ki++]=(double)i;
hud[ki++]=a1->vx;
hud[ki++]=a1->vy;
hud[ki++]=a1->vz;
// printf("rank=%i i=%i ax=%le ay=%le az=%le vx=%le \n",rank,i,a1->ax,a1->ay,a1->az,a1->vx);
}
if(ki == 1) ki=0;
if(kj == 1) kj=0;
hud[0]=(double)ki;
hvd[0]=(double)kj;
//printf("rankmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm=%i kj=%i ki=%i \n",rank,kj,ki);
if(numprocs > 1){
MPI_Sendrecv(&hud[0],ki,MPI_DOUBLE,left,5,&had[0],5000,MPI_DOUBLE,right,5,comm1d,&status);
MPI_Sendrecv(&hvd[0],kj,MPI_DOUBLE,right,6,&hbd[0],5000,MPI_DOUBLE,left,6,comm1d,&status);
}
ki = (int)had[0];
kj = 1;
////printf("rankrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr=%i ki=%i \n",rank,ki);
while( kj < ki){
i=(int)had[kj++];
a1=at+i;
if( a1->typ != 1)
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRR2");
printf(" rank=%i i=%i a1->typ=%i a1->xreal=%le xsl=%i xer=%i \n",rank,i,a1->typ,a1->xreal,xsl,xer);
exit(1);
}
a1->vx=had[kj++];
a1->vy=had[kj++];
a1->vz=had[kj++];
//printf("rank=%i i=%i ax=%le ay=%le az=%le vx=%le \n",rank,i,a1->ax,a1->ay,a1->az,a1->vx);
// printf(" %i\t%i\t%i\t%le\t%i\t%i\n",ki,rank,a->typ,a->xreal,xsl,xer);
}
ki = (int)hbd[0];
kj = 1;
//printf("rankffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff=%i ki=%i \n",rank,ki);
while( kj < ki){
i=(int)hbd[kj++];
a1=at+i;
if( a1->typ != 1)
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRR3");
printf(" rank=%i i=%i a1->typ=%i a1->xreal=%le xsl=%i xer=%i \n",rank,i,a1->typ,a1->xreal,xsl,xer);
exit(1);
}
a1->vx=hbd[kj++];
a1->vy=hbd[kj++];
a1->vz=hbd[kj++];
// printf("rank=%i i=%i ax=%le ay=%le az=%le vx=%le \n",rank,i,a1->ax,a1->ay,a1->az,a1->vx);
}
for(i=0;i<Naf;i++)
{
a1=at+i;
if((a1->typ == 1) && (a1->xreal >= xsl) && (a1->xreal < xer))
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRR4");
printf(" rank=%i i=%i a1->typ=%i a1->xreal=%le xsl=%i xer=%i \n",rank,i,a1->typ,a1->xreal,xsl,xer);
exit(1);
}
if((a1->typ == 0) && ((a1->xreal < xsl) || (a1->xreal >= xer)))
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRR5");
printf(" rank=%i i=%i a1->typ=%i a1->xreal=%le xsl=%i xer=%i \n",rank,i,a1->typ,a1->xreal,xsl,xer);
exit(1);
}
}
for(i=0;i<length;i++)
{
a1=at+i;
if( a1->typ == 2)continue;
if(((a1->xreal) < (xer+1.4*(xer-xsl))) && ((a1->xreal) > (xsl-1.4*(xer-xsl))))
{
k=i;
j=0;
dax=day=daz=0;
c=cnt+i;
do
{
a1=(at+(k+(j*length)));
dax+=a1->xreal;
day+=a1->yreal;
daz+=a1->zreal;
//c->bspring[j]=a1->n;
a1->plane=k;
j++;
}
while(j<npoints);
dax*=inv_npoints;
day*=inv_npoints;
daz*=inv_npoints;
c->xreal=dax;
c->yreal=day;
c->zreal=daz;
c->radius=0;
c->radius2=0;
}
}
//exit(1);
for(i=1;i<length;i++)
{
a1=at+i;
if( a1->typ == 2)continue;
if((((cnt+i)->xreal) < (xer+(xer-xsl))) && (((cnt+i)->xreal) > (xsl-(xer-xsl))))
{
dax=(cnt+i-1)->xreal-(cnt+i)->xreal;
day=(cnt+i-1)->yreal-(cnt+i)->yreal;
daz=(cnt+i-1)->zreal-(cnt+i)->zreal;
inv_d=1./sqrt(dax*dax+day*day+daz*daz);
(cnt+i-1)->dirx=dax*inv_d;
(cnt+i-1)->diry=day*inv_d;
(cnt+i-1)->dirz=daz*inv_d;
if(i==length-1)
{
inv_d=1./sqrt(dax*dax+day*day+daz*daz);
(cnt+i)->dirx=dax*inv_d;
(cnt+i)->diry=day*inv_d;
(cnt+i)->dirz=daz*inv_d;
}
}
}
for(i=0;i<Naf;i++)
{
a1=at+i;
if(a1->xreal < (xsl-4)) a1->typ=2;
if(a1->xreal > (xer+4)) a1->typ=2;
}
for(i=0;i<(length*npoints);i++)
{
a1=at+i;
if(a1->typ == 2)continue;
c=cnt+a1->plane;
dax=c->xreal-a1->xreal;
day=c->yreal-a1->yreal;
daz=c->zreal-a1->zreal;
d=sqrt(dax*dax+day*day+daz*daz);
if(d>0)
{
inv_d=1./d;
a1->nx=-dax*inv_d;
a1->ny=-day*inv_d;
a1->nz=-daz*inv_d;
if(tsim==0)
{
a1->nx0=-dax*inv_d;
a1->ny0=-day*inv_d;
a1->nz0=-daz*inv_d;
//printf("%le\t%le\t%le\n",a1->xreal+a1->nx0,a1->yreal+a1->ny0,a1->zreal+a1->nz0);
}
if(d>c->radius)
{
c->radius2=dax*dax+day*day+daz*daz;
c->radius=d;
}
}
}
if(svent==1)
{
for(i=0;i<Naf;i++)
{
a1=at+i;
if(a1->typ == 2)continue;
a1->cx=a1->xreal;
a1->cy=a1->yreal;
a1->cz=a1->zreal;
//printf("%le\t%le\t%le\n",a1->cx,a1->cy,a1->cz);
}
for(i=0;i<length;i++)
{
c=cnt+i;
c->cx=c->xreal;
c->cy=c->yreal;
c->cz=c->zreal;
//printf("%le\t%le\t%le\n",a1->cx,a1->cy,a1->cz);
}
for(i=0;i<Naf;i++)
{
a1=at+i;
if(a1->xreal < (xsl-3)) a1->typ=2;
if(a1->xreal > (xer+3)) a1->typ=2;
}
}
//exit(1);
}
void rbcpos()
{
unsigned long int i,ii,j,n,ki,kj,il,idd;
int k;
double xl,yl,zl,rpc,dcp;
struct redcell *rc;
struct aftr *a;
///////////
static double *pu; ///////////////
static double *pra; //////////////
static int n1st = 0;
if(n1st == 0){
pu = (double *) malloc(sizeof(double)*(5000));////////////////
if(pu == NULL){
printf("\n err_message(1)");/////////////////////
exit(1);
}
pra = (double *) malloc(sizeof(double)*(5000));/////////////////
if(pra == NULL) exit(1);////////////parallel///////////
n1st++;
}
//////////
Nrbcp=0;
if(OBSTYP==0)
{
printf(" OBSTYP=%i\n",OBSTYP);
j=0;
i=0;
k=0;
n=0;
Nrbct=1650;
Rrbc=1.0;
Mrbc=100*nprcell*4*PI*Rrbc*Rrbc*Rrbc/3;
Irbc=2*Mrbc*Rrbc*Rrbc/5;
if(rank==0)
{
do
{
// xl=(double)xsl+drand48()*(double)(L/numprocs);
xl=drand48()*(double)L;
yl= drand48()*(double)Lr1;
zl= drand48()*(double)Lr1;
rpc=sqrt((yl-ylc)*(yl-ylc)+(zl-zlc)*(zl-zlc));
if( (rpc >= (Lr0-Rrbc)) || (xl <= 2*Rrbc) || (xl >= ((double)L-2*Rrbc)) ) continue; //
// if( (rank< 5 && rpc < (1.4+Rrbc)) || (rpc >= (Lr0-Rrbc)) || (xl <= (xsl-Rrbc)) || (xl >= (xer+Rrbc)) ) continue;
for(i=0;i<j;i++)
{
rc=rbc+i;
dcp=sqrt((xl-rc->xcm)*(xl-rc->xcm)+(yl-rc->ycm)*(yl-rc->ycm)+(zl-rc->zcm)*(zl-rc->zcm));
if(dcp <= 2.05*Rrbc) i=Nrbct+10;
}
for(ii=0;ii<Naf;ii++)
{
a=at+ii;
dcp=sqrt((xl-a->xreal)*(xl-a->xreal)+(yl-a->yreal)*(yl-a->yreal)+(zl-a->zreal)*(zl-a->zreal));
if(dcp <= 1.1*Rrbc) ii=Naf+10;
}
if(i < Nrbct+2 && ii< Naf+2)
{
rc=rbc+j;
rc->id=j;
rc->xcm=xl;
rc->ycm=yl;
rc->zcm=zl;
/* rc->xcm=15.0+j*10;
if(j < 5)
{
rc->ycm=11.0;
rc->zcm=10.0;
}
else
{
rc->ycm=8.0;
rc->zcm=8.0;
}*/
j++;
printf(" %i\t%i\t%i\t%le\t%le\t%le\n",rank,j,i,rc->xcm,rc->ycm,rc->zcm);
}
}
while(j< Nrbct);
for(i=0;i<Nrbct;i++)
{
rc=rbc+i;
printf(" %i\t%i\t%i\t%le\t%le\t%le\n",rank,j,i,rc->xcm,rc->ycm,rc->zcm);
}
}
for(k=1;k<numprocs;k++)
{
pu[0] = 0.0;
ki=1;
if(rank==0)
{
for(i=0;i<j;i++)
{
rc=rbc+i;
xl=rc->xcm;
if(xl >= (k*(L/numprocs)-Rrbc-0.5) && (xl <= ((k+1)*(L/numprocs)+Rrbc+0.5)) )
{
pu[ki++]=(double)rc->id;
pu[ki++]=rc->xcm;
pu[ki++]=rc->ycm;
pu[ki++]=rc->zcm;
}
}
if(ki == 1) ki=0;
pu[0]=(double)ki;
MPI_Send(&pu[0],ki,MPI_DOUBLE,k,k,MPI_COMM_WORLD);
}
if(rank==k)
{
MPI_Recv(&pra[0],5000,MPI_DOUBLE,0,rank,MPI_COMM_WORLD,&status);
ki = (int)pra[0];
kj = 1;
il=0;
while( kj < ki)
{
rc=rbc+il;
rc->id=(int)pra[kj++];
rc->xcm=pra[kj++];
rc->ycm=pra[kj++];
rc->zcm=pra[kj++];
il++;
}
Nrbcp=il;
}
}
if(rank==0)
{
il=0;
for(i=0;i<j;i++)
{
rc=rbc+i;
idd=rc->id;
xl=rc->xcm;
yl=rc->ycm;
zl=rc->zcm;
if( xl <= ((L/numprocs)+Rrbc+0.5) )
{
rc=rbc+il;
rc->id=idd;
rc->xcm=xl;
rc->ycm=yl;
rc->zcm=zl;
il++;
}
}
Nrbcp=il;
}
for(i=0;i<Nrbcp;i++)
{
rc=rbc+i;
rc->ip=0;
xl=rc->xcm;
if( (xl <= (xsl+Rrbc+0.5)) || (xl >= (xer-Rrbc-0.5)) )
{
rc->ip=1;
}
}
for(i=0;i<Nrbcp;i++)
{
rc=rbc+i;
rc->Vxcm=0.0;
rc->Vycm=0.0;
rc->Vzcm=0.0;
rc->wx=0.0;
rc->wy=0.0;
rc->wz=0.0;
rc->ex=1.0;
rc->ey=0.0;
rc->ez=0.0;
rc->Momx=0.0;
rc->Momy=0.0;
rc->Momz=0.0;
rc->Amomx=0.0;
rc->Amomy=0.0;
rc->Amomz=0.0;
printf(" %i\t%i\t%i\t%i\t%i\t%le\t%le\t%le\n",rank,Nrbcp,i,rc->id,rc->ip,rc->xcm,rc->ycm,rc->zcm);
}
}
}
void solventposvel1()
{
unsigned long int i,j,k,fsolvent,n;
double xl,yl,zl,v2,Vx,Vy,Vz,temp,inv_temp,rpw,rpc;
struct sol *p;
struct wallp *pw;
struct redcell *rc;
srand48(rank*100*Nsol);
i=0;
k=0;
j=0;
n=0;
do
{
xl=(double)xsl+drand48()*(double)(L/numprocs);
yl=drand48()*(double)Lr1;
zl=drand48()*(double)Lr1;
rpw=sqrt((yl-ylc)*(yl-ylc)+(zl-zlc)*(zl-zlc));
if(rpw > Lr0)
{
pw=ww+k;
pw->x=xl;
pw->y=yl;
pw->z=zl;
pw->n=k;
pw->vx=gaussianrandom();
pw->vy=gaussianrandom();
pw->vz=gaussianrandom();
k++;
}
else
{
rpc=100.0;
for(i=0;i<Nrbcp;i++)
{
rc=rbc+i;
if(rc->ip == 2)continue;
rpc=sqrt((xl-rc->xcm)*(xl-rc->xcm)+(yl-rc->ycm)*(yl-rc->ycm)+(zl-rc->zcm)*(zl-rc->zcm));
if(rpc<Rrbc)
{
//if(rank==0) printf(" %i\t%i\t%i\t%le\t%le\t%le\t%le\n",rank,Nrbcp,i,rc->xcm,rc->ycm,rc->zcm,rpc);
i=50000;
}
}
if(rpc < Rrbc)
{
pw=ww+k;
pw->x=xl;
pw->y=yl;
pw->z=zl;
pw->n=k;
pw->vx=gaussianrandom();
pw->vy=gaussianrandom();
pw->vz=gaussianrandom();
k++;
n++;
}
else
{
p=ss+j;
p->x=xl;
p->y=yl;
p->z=zl;
p->n=j;
j++;
}
}
}
while((j+k)<Nsolp);
Nsolp=j;
Nsolw=k;
printf("rank=%i nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnN=%i k=%i j=%i \n",rank,n,k,j);
v2=Vx=Vy=Vz=0;
for(i=0;i<Nsolp;i++)
{
p=ss+i;
p->vx=gaussianrandom();
p->vy=gaussianrandom();
p->vz=gaussianrandom();
if(i < 100) printf("rank=%i i=%i p->vx=%le \n",rank,i,p->vx);
Vx+=p->vx;
Vy+=p->vy;
Vz+=p->vz;
v2+=p->vx*p->vx+p->vy*p->vy+p->vz*p->vz;
}
printf("rank=%i Nsol=%i Nsolp=%i sol inside af solTEMP=%le\n",rank,Nsol,Nsolp,1./3.*v2/Nsolp);
//}
}
void wallpos()
{
unsigned long int i,j,k,n;
double xl,yl,zl,rpw,rpc;
double Vxrot,Vyrot,Vzrot,vsx,vsy,vsz;
struct wallp *pw;
struct redcell *rc;
k=0;
j=0;
n=0;
do
{
xl=(double)xsl-0.5+drand48()*((double)(L/numprocs)+1.0);
yl=drand48()*(double)Lr1;
zl=drand48()*(double)Lr1;
rpw=sqrt((yl-ylc)*(yl-ylc)+(zl-zlc)*(zl-zlc));
if(rpw > Lr0)
{
pw=ww+k;
pw->x=xl;
pw->y=yl;
pw->z=zl;
pw->n=k;
pw->vx=gaussianrandom();
pw->vy=gaussianrandom();
pw->vz=gaussianrandom();
k++;
}
else
{
rpc=100.0;
for(i=0;i<Nrbcp;i++)
{
rc=rbc+i;
if(rc->ip == 2)continue;
rpc=sqrt((xl-rc->xcm)*(xl-rc->xcm)+(yl-rc->ycm)*(yl-rc->ycm)+(zl-rc->zcm)*(zl-rc->zcm));
if(rpc<Rrbc)
{
//if(rank==0) printf(" %i\t%i\t%i\t%le\t%le\t%le\t%le\n",rank,Nrbcp,i,rc->xcm,rc->ycm,rc->zcm,rpc);
//////////
pw=ww+k;
pw->x=xl;
pw->y=yl;
pw->z=zl;
pw->n=k;
pw->vx=gaussianrandom();
pw->vy=gaussianrandom();
pw->vz=gaussianrandom();
if(OBSM==1)
{
Vxrot=crossx(rc->wx,rc->wy,rc->wz,(xl-rc->xcm),(yl-rc->ycm),(zl-rc->zcm));
Vyrot=crossy(rc->wx,rc->wy,rc->wz,(xl-rc->xcm),(yl-rc->ycm),(zl-rc->zcm));
Vzrot=crossz(rc->wx,rc->wy,rc->wz,(xl-rc->xcm),(yl-rc->ycm),(zl-rc->zcm));
pw->vx+=rc->Vxcm+Vxrot;
pw->vy+=rc->Vycm+Vyrot;
pw->vz+=rc->Vzcm+Vzrot;
}
k++;
n++;
i=50000;
}
}
}
j++;
}
while(j < ((Nsol/numprocs)+(Nsol/L)));
Nsolw=k;
//printf("rank=%i nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnN=%i k=%i j=%i \n",rank,n,k,j);
}
void streaming(double stime)
{
unsigned long int i,ic,j,out,n,ki=0,kj=0,k2=0,k3=0;
double xl,yl,zl,rnp,v2,x1,y1,z1,st,st1,rpc;
struct sol *srd;
struct redcell *rc;
/////////////////////////////
static double *hup;
static double *hvp;
static double *hap;
static double *hbp;
static int *savid;
static int n2st = 0;
if(n2st == 0){
hup = (double *) malloc(sizeof(double)*(600000));
if(hup == NULL){
printf("\n err_message(2)");
exit(1);
}
hvp = (double *) malloc(sizeof(double)*(600000));
if(hvp == NULL) exit(1);/////////
hap = (double *) malloc(sizeof(double)*(600000));
if(hap == NULL) exit(1);/////////
hbp = (double *) malloc(sizeof(double)*(600000));
if(hbp == NULL) exit(1);/////////
savid = (int *) malloc(sizeof(int)*(200000));
if(savid == NULL) exit(1);//////////
savid[0]=1;
n2st++;
}
// MPI_Status status;
hup[ki++] = 0.0;
hvp[kj++] = 0.0;
hap[0] = 0.0;
hbp[0] = 0.0;
k2=savid[0];
/////////////////////////////////////////
//printf("ranka=%i k2=%i savid[0]=%i kj=%i ki=%i \n",rank,k2,savid[0],kj,ki);
v2=0;
for(i=0;i<Nsolp;i++)
{
srd=ss+i;
///////////////////////////////
xl=srd->x;
if(xl < -1000)
continue;
//////////////////////////////////////////
srd->in=0;
out=j=0;
st1=st=stime;
do
{
xl=srd->x+srd->vx*st;
yl=srd->y+srd->vy*st;
zl=srd->z+srd->vz*st;
///printf("rank=%i n=%i x=%le y=%le z=%le xc=%le yc=%le zc=%le r2=%le radius2=%le \n",rank,srd->n,xl,yl,zl,(cnt+length-1)->cx,(cnt+length-1)->cy,(cnt+length-1)->cz,r2,(cnt+length-1)->radius2);
if(numprocs == 1){
if(xl<0)
xl+=(double)L;
if(xl>=(double)L)
xl-=(double)L;
}
/*
if(yl<0)
yl+=(double)L;
if(yl>=(double)L)
yl-=(double)L;
if(zl<0)
zl+=(double)L;
if(zl>=(double)L)
zl-=(double)L;
*/
srd->x=xl;
srd->y=yl;
srd->z=zl;
if(j<10)
{
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
rpc=sqrt((xl-rc->xcm)*(xl-rc->xcm)+(yl-rc->ycm)*(yl-rc->ycm)+(zl-rc->zcm)*(zl-rc->zcm));
if(rpc<Rrbc)
{
// if(rank==0) printf(" %i\t%i\t%i\t%le\t%le\t%le\t%le\n",rank,Nrbcp,i,rc->xcm,rc->ycm,rc->zcm,rpc);
bouncebackrbc(0,i,ic,st);
ic=50000;
}
}
}
////////////////////
//insideaftr8(srd->n);
// yl=srd->y;
// zl=srd->z;
rnp=sqrt((yl-ylc)*(yl-ylc)+(zl-zlc)*(zl-zlc));
if(rnp > Lr0)
{
bounceback(0,i,st);
}
j++;
insideaftr8(srd->n);
/* if(srd->in != 0 && tsim == 0)
{
srd->x=-2000.0;
srd->in==0;
printf("ttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE rank=%i i=%i \n",rank,i);
}
*/
st1*=0.5;
st=st1;
if(srd->in==0 || j>20)
out=1;
}
while(out==0);
if((srd->x > xer+2) || (srd->x < xsl-2))
{
printf("EEEEEEEEEEEEEEEEEEEEEEEEEERRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRROOOOOOORRRRRRRR rank=%i i=%i \n",rank,i);
srd->x=-2000.0;
continue;
}
yl=srd->y;
zl=srd->z;
rnp=sqrt((yl-ylc)*(yl-ylc)+(zl-zlc)*(zl-zlc));
if(rnp > (Lr0+0.2))
{
printf("EEEEEEEEEEEEEEERRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRROOOOOOORRRRRRRR rank=%i i=%i \n",rank,i);
srd->x=-2000.0;
continue;
}
//if(j > 1) printf("rank=%i j=%i \n",rank,j);
/////////////////////////////////////////
//k2=1;//savid[0];
if(srd->x >= xer){
hvp[kj++]=srd->x;
hvp[kj++]=srd->y;
hvp[kj++]=srd->z;
hvp[kj++]= srd->vx;
hvp[kj++]= srd->vy;
hvp[kj++]= srd->vz;
savid[k2++]=srd->n;
srd->x=-2000.0;
}
else if(srd->x < xsl){
hup[ki++]=srd->x;
hup[ki++]=srd->y;
hup[ki++]=srd->z;
hup[ki++]= srd->vx;
hup[ki++]= srd->vy;
hup[ki++]= srd->vz;
savid[k2++]=srd->n;
srd->x=-2000.0;
}
/////////////////////////////////////////////////////
}
savid[0]=k2;
if(ki == 1) ki=0;
if(kj == 1) kj=0;
hup[0]=(double)ki;
hvp[0]=(double)kj;
if(numprocs > 1){
MPI_Sendrecv(&hup[0],ki,MPI_DOUBLE,left,3,&hap[0],600000,MPI_DOUBLE,right,3,comm1d,&status);
MPI_Sendrecv(&hvp[0],kj,MPI_DOUBLE,right,4,&hbp[0],600000,MPI_DOUBLE,left,4,comm1d,&status);
}
ki = (int)hap[0];
kj = 1;
k3=1;
while( kj < ki){
if(k3 < savid[0]){
j = savid[k3++];
}
else{
j = Nsolp ;
Nsolp++;
printf("rank=%i Nsolp=%i savid[0]=%i \n",rank,Nsolp,savid[0]);
}
if(j < 0)
{
printf("EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEERRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRROOOOOOORRRRRRRR rank=%i j=%i \n",rank,j);
exit(1);
}
srd=ss+j;
srd->n=j;
srd->x=hap[kj++];
if(rank==numprocs-1) srd->x+=(double)L;
srd->y=hap[kj++];
srd->z=hap[kj++];
srd->vx=hap[kj++];
srd->vy=hap[kj++];
srd->vz=hap[kj++];
}
ki = (int)hbp[0];
kj = 1;
while( kj < ki){
if(k3 < savid[0] ){
j = savid[k3++];
}
else{
j = Nsolp ;
Nsolp++;
printf("rank1111111=%i Nsolp=%i savid[0]=%i \n",rank,Nsolp,savid[0]);
}
if(j < 0)
{
printf("2222222222EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEERRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRROOOOOOORRRRRRRR rank=%i j=%i \n",rank,j);
exit(1);
}
srd=ss+j;
srd->n=j;
srd->x=hbp[kj++];
if(rank == 0) srd->x-=(double)L;
srd->y=hbp[kj++];
srd->z=hbp[kj++];
srd->vx=hbp[kj++];
srd->vy=hbp[kj++];
srd->vz=hbp[kj++];
}
j=1;
for(i=k3;i< savid[0];i++)
{
savid[j++]=savid[i];
}
savid[0]=j;
for(i=0;i<Nsolp;i++)
{
srd=ss+i;
/////////////////////////parallel///////////
xl=srd->x;
if(xl < -1000)
continue;
if(xl < 0 || srd->y < 0 || srd->z < 0 || xl >= L || srd->y >= L || srd->z >= L )
{
printf("EEEEEEEEEEEEEEEEEEEEERRRRRRRRRRRRRRRRRRRRRRRROOOOOOORRRRRRRR rank=%i xl=%le \n",rank,xl);
exit(1);
}
if(xl < xsl || xl >= xer)
{
printf("EEEEEEEEEEEEEEEEEEEEERRRRRRRRRRRRRRRRRRRRRRRROOOOOOORRRRRRRR11111111 rank=%i i=%i xl=%le \n",rank,i,xl);
exit(1);
}
}
}
void bounceback(int k ,unsigned long int n, double st)
{
double x1,y1,z1,ey,ez,r,d,r2;
struct sol *srd;
struct aftr *a;
if(k == 0)
{
srd=ss+n;
x1=srd->x;
y1=srd->y;
z1=srd->z;
if(st < 0.001)
{
///printf("EEEEERRRRRRROOOOOOORRRRRRRRRRRRRRdt rank=%ld x=%le y=%le z=%le tsim=%ld n=%ld st=%le \n",rank,x1,y1,z1,tsim,n,st);
}
x1=x1-srd->vx*st*0.5;
y1=y1-srd->vy*st*0.5;
z1=z1-srd->vz*st*0.5;
r=sqrt((y1-ylc)*(y1-ylc)+(z1-zlc)*(z1-zlc));
d=Lr0-r;
ey=(y1-ylc)/r;
ez=(z1-zlc)/r;
y1=y1+d*ey;
z1=z1+d*ez;
srd->vx=-srd->vx;
srd->vy=-srd->vy;
srd->vz=-srd->vz;
srd->x=x1+srd->vx*st*0.5;
srd->y=y1+srd->vy*st*0.5;
srd->z=z1+srd->vz*st*0.5;
}
else if(k ==1)
{
if(st > 0.001)
{
printf("EEEEERRRRRRRRRRRRRRRRRRRRRRRRROOOOOOOOOOOOOOOOOOOOOOOOOOOOOORRRRRRRRRRRRRRRRdt tsim=%ld n=%ld st=%le \n",tsim,n,st);
exit(1);
}
a=at+n;
x1=a->xreal;
y1=a->yreal;
z1=a->zreal;
x1=x1-a->vx*st*0.5;
y1=y1-a->vy*st*0.5;
z1=z1-a->vz*st*0.5;
r=sqrt((y1-ylc)*(y1-ylc)+(z1-zlc)*(z1-zlc));
d=Lr0-r;
ey=(y1-ylc)/r;
ez=(z1-zlc)/r;
y1=y1+d*ey;
z1=z1+d*ez;
a->vx=a->vx;
a->vy=-a->vy;
a->vz=-a->vz;
a->xreal=x1+a->vx*st*0.5;
a->yreal=y1+a->vy*st*0.5;
a->zreal=z1+a->vz*st*0.5;
}
}
void bouncebackrbc(int k ,unsigned long int n,int ic, double st)
{
double x1,y1,z1,ex,ey,ez,r,d,r2;
double Vxrot,Vyrot,Vzrot,vsx,vsy,vsz;
struct sol *srd;
struct aftr *a;
struct redcell *rc;
if(k == 0)
{
srd=ss+n;
x1=srd->x;
y1=srd->y;
z1=srd->z;
x1=x1-srd->vx*st*0.5;
y1=y1-srd->vy*st*0.5;
z1=z1-srd->vz*st*0.5;
rc=rbc+ic;
if(OBSTYP==0)
{
r=sqrt((x1-rc->xcm)*(x1-rc->xcm)+(y1-rc->ycm)*(y1-rc->ycm)+(z1-rc->zcm)*(z1-rc->zcm));
d=Rrbc-r;
ex=(x1-rc->xcm)/r;
ey=(y1-rc->ycm)/r;
ez=(z1-rc->zcm)/r;
x1=x1+d*ex;
y1=y1+d*ey;
z1=z1+d*ez;
if(OBSM==0)
{
srd->vx=-srd->vx;
srd->vy=-srd->vy;
srd->vz=-srd->vz;
}
else if(OBSM==1)
{
Vxrot=crossx(rc->wx,rc->wy,rc->wz,(x1-rc->xcm),(y1-rc->ycm),(z1-rc->zcm));
Vyrot=crossy(rc->wx,rc->wy,rc->wz,(x1-rc->xcm),(y1-rc->ycm),(z1-rc->zcm));
Vzrot=crossz(rc->wx,rc->wy,rc->wz,(x1-rc->xcm),(y1-rc->ycm),(z1-rc->zcm));
vsx=srd->vx;
vsy=srd->vy;
vsz=srd->vz;
srd->vx=-srd->vx+rc->Vxcm+Vxrot;
srd->vy=-srd->vy+rc->Vycm+Vyrot;
srd->vz=-srd->vz+rc->Vzcm+Vzrot;
rc->Momx+=srd->vx-vsx;
rc->Momy+=srd->vy-vsy;
rc->Momz+=srd->vz-vsz;
rc->Amomx+=crossx((x1-rc->xcm),(y1-rc->ycm),(z1-rc->zcm),(srd->vx-vsx),(srd->vy-vsy),(srd->vz-vsz));
rc->Amomy+=crossy((x1-rc->xcm),(y1-rc->ycm),(z1-rc->zcm),(srd->vx-vsx),(srd->vy-vsy),(srd->vz-vsz));
rc->Amomz+=crossz((x1-rc->xcm),(y1-rc->ycm),(z1-rc->zcm),(srd->vx-vsx),(srd->vy-vsy),(srd->vz-vsz));
}
srd->x=x1+srd->vx*st*0.5;
srd->y=y1+srd->vy*st*0.5;
srd->z=z1+srd->vz*st*0.5;
//r2=sqrt((srd->x-rc->xcm)*(srd->x-rc->xcm)+(srd->y-rc->ycm)*(srd->y-rc->ycm)+(srd->z-rc->zcm)*(srd->z-rc->zcm));
//if(rank==0) printf(" %i\t%i\t%i\t%le\t%le\t%le\t%le\n",rank,Nrbcp,n,rc->xcm,rc->ycm,rc->zcm,r2);
}
}
else if(k == 1)
{
a=at+n;
x1=a->xreal;
y1=a->yreal;
z1=a->zreal;
x1=x1-a->vx*st*0.5;
y1=y1-a->vy*st*0.5;
z1=z1-a->vz*st*0.5;
rc=rbc+ic;
if(OBSTYP==0)
{
r=sqrt((x1-rc->xcm)*(x1-rc->xcm)+(y1-rc->ycm)*(y1-rc->ycm)+(z1-rc->zcm)*(z1-rc->zcm));
d=Rrbc-r;
ex=(x1-rc->xcm)/r;
ey=(y1-rc->ycm)/r;
ez=(z1-rc->zcm)/r;
x1=x1+d*ex;
y1=y1+d*ey;
z1=z1+d*ez;
if(OBSM==0)
{
a->vx=a->vx;
a->vy=-a->vy;
a->vz=-a->vz;
}
else if(OBSM==1)
{
Vxrot=crossx(rc->wx,rc->wy,rc->wz,(x1-rc->xcm),(y1-rc->ycm),(z1-rc->zcm));
Vyrot=crossy(rc->wx,rc->wy,rc->wz,(x1-rc->xcm),(y1-rc->ycm),(z1-rc->zcm));
Vzrot=crossz(rc->wx,rc->wy,rc->wz,(x1-rc->xcm),(y1-rc->ycm),(z1-rc->zcm));
vsx=a->vx;
vsy=a->vy;
vsz=a->vz;
a->vx=a->vx+rc->Vxcm+Vxrot;
a->vy=-a->vy+rc->Vycm+Vyrot;
a->vz=-a->vz+rc->Vzcm+Vzrot;
rc->Momx+=a->vx-vsx;
rc->Momy+=a->vy-vsy;
rc->Momz+=a->vz-vsz;
rc->Amomx+=crossx((x1-rc->xcm),(y1-rc->ycm),(z1-rc->zcm),(a->vx-vsx),(a->vy-vsy),(a->vz-vsz));
rc->Amomy+=crossy((x1-rc->xcm),(y1-rc->ycm),(z1-rc->zcm),(a->vx-vsx),(a->vy-vsy),(a->vz-vsz));
rc->Amomz+=crossz((x1-rc->xcm),(y1-rc->ycm),(z1-rc->zcm),(a->vx-vsx),(a->vy-vsy),(a->vz-vsz));
}
a->xreal=x1+a->vx*st*0.5;
a->yreal=y1+a->vy*st*0.5;
a->zreal=z1+a->vz*st*0.5;
// r2=sqrt((a->xreal-rc->xcm)*(a->xreal-rc->xcm)+(a->yreal-rc->ycm)*(a->yreal-rc->ycm)+(a->zreal-rc->zcm)*(a->zreal-rc->zcm));
// printf(" %i\t%i\t%i\t%le\t%le\t%le\t%le\n",rank,Nrbcp,n,rc->xcm,rc->ycm,rc->zcm,r2);
}
}
}
void insideaftr8(unsigned long int n)
{
unsigned long int i,w1,w2;
double x1,y1,z1,d2,r2,inv_d;
double nx,ny,nz,tx1,ty1,tz1,dsmall,ds1;
double r,rm1,rp1,dirvec,trvx,trvy,trvz,scalar;
struct sol *srd;
struct cnline *c;
struct aftr *a,*a1;
srd=ss+n;
srd->in=0;
w1=5000;
dsmall=30.;
for(i=0;i<length;i++)
{
a1=at+i;
if(a1->typ == 2 )continue;
d2=dist2(a1->cx,srd->x,a1->cy,srd->y,a1->cz,srd->z);
if(d2<dsmall)
{
x1=srd->x;
y1=srd->y;
z1=srd->z;
w1=i;
dsmall=d2;
}
}
w2=Naf+2;
if(w1<length)
{
ds1=30.;
for(i=0;i<Naf;i++)
{
a1=at+i;
if(a1->typ == 2)continue;
d2=dist2(a1->cx,x1,a1->cy,y1,a1->cz,z1);
if(d2<ds1 && d2<6.)
{
ds1=d2;
w2=i;
a=at+w2;
scalar=a->nx*(a->cx-x1)+a->ny*(a->cy-y1)+a->nz*(a->cz-z1);
/// scalar1=a->nx*(srd->vx-a->vx)+a->ny*(srd->vy-a->vy)+a->nz*(srd->vz-a->vz);
}
}
}
//if((w2>24 && w2<30) || (w2>54 && w2<60) || (w2>84 && w2<90) || (w2>114 && w2<120) || (w2>144 && w2<150) || (w2>174 && w2<180) || (w2>204 && w2<210)) w2=2000;
//if((w2>234 && w2<240) || (w2>264 && w2<270) || (w2>294 && w2<300)) w2=2000;
if(w2<Naf && scalar>0)
{
a=at+w2;
//printf("%le\t%le\t%le\t",x1,y1,z1);
if(a->plane!=0 && a->plane!=length-1)
{
trvx=(cnt+a->plane)->dirx;
trvy=(cnt+a->plane)->diry;
trvz=(cnt+a->plane)->dirz;
dirvec=trvx*(a->cx-x1)+trvy*(a->cy-y1)+trvz*(a->cz-z1);
if(dirvec<0)
{
r=(cnt+a->plane)->radius;
rm1=(cnt+(at+w2-1)->plane)->radius;
if(r-rm1<0)
{
nx=(at+w2-1)->nx;
ny=(at+w2-1)->ny;
nz=(at+w2-1)->nz;
scalar=nx*((at+w2-1)->cx-x1)+ny*((at+w2-1)->cy-y1)+nz*((at+w2-1)->cz-z1);
srd->x=x1+nx*scalar;
srd->y=y1+ny*scalar;
srd->z=z1+nz*scalar;
velupdate2(a->n,srd->n);
return;
}
else
{
nx=a->nx;
ny=a->ny;
nz=a->nz;
scalar=a->nx*(a->cx-x1)+a->ny*(a->cy-y1)+a->nz*(a->cz-z1);
srd->x=x1+nx*scalar;
srd->y=y1+ny*scalar;
srd->z=z1+nz*scalar;
velupdate2(a->n,srd->n);
return;
}
}
else
{
r=(cnt+a->plane)->radius;
rp1=(cnt+(at+w2+1)->plane)->radius;
if(r-rp1<0)
{
nx=(at+w2+1)->nx;
ny=(at+w2+1)->ny;
nz=(at+w2+1)->nz;
scalar=nx*((at+w2+1)->cx-x1)+ny*((at+w2+1)->cy-y1)+nz*((at+w2+1)->cz-z1);
srd->x=x1+nx*scalar;
srd->y=y1+ny*scalar;
srd->z=z1+nz*scalar;
velupdate2(a->n,srd->n);
return;
}
else
{
nx=a->nx;
ny=a->ny;
nz=a->nz;
scalar=a->nx*(a->cx-x1)+a->ny*(a->cy-y1)+a->nz*(a->cz-z1);
srd->x=x1+nx*scalar;
srd->y=y1+ny*scalar;
srd->z=z1+nz*scalar;
velupdate2(a->n,srd->n);
return;
}
}
}
if(a->plane==0)
{
trvx=(cnt+a->plane)->dirx;
trvy=(cnt+a->plane)->diry;
trvz=(cnt+a->plane)->dirz;
dirvec=trvx*(a->cx-x1)+trvy*(a->cy-y1)+trvz*(a->cz-z1);
if(dirvec>0)
{
nx=(at+w2+1)->nx;
ny=(at+w2+1)->ny;
nz=(at+w2+1)->nz;
scalar=nx*((at+w2+1)->cx-x1)+ny*((at+w2+1)->cy-y1)+nz*((at+w2+1)->cz-z1);
srd->x=x1+nx*scalar;
srd->y=y1+ny*scalar;
srd->z=z1+nz*scalar;
velupdate2(a->n,srd->n);
return;
}
else
{
r2=dist2((cnt+0)->cx,x1,(cnt+0)->cy,y1,(cnt+0)->cz,z1);
if(r2<(cnt+0)->radius2)
{
nx=x1-(cnt+0)->cx;
ny=y1-(cnt+0)->cy;
nz=z1-(cnt+0)->cz;
inv_d=1./sqrt(nx*nx+ny*ny+nz*nz);
nx*=inv_d;
ny*=inv_d;
nz*=inv_d;
tx1=crossx(nx,ny,nz,((cnt+0)->cx-a->cx),((cnt+0)->cy-a->cy),((cnt+0)->cz-a->cz));
ty1=crossy(nx,ny,nz,((cnt+0)->cx-a->cx),((cnt+0)->cy-a->cy),((cnt+0)->cz-a->cz));
tz1=crossz(nx,ny,nz,((cnt+0)->cx-a->cx),((cnt+0)->cy-a->cy),((cnt+0)->cz-a->cz));
inv_d=1./sqrt(tx1*tx1+ty1*ty1+tz1*tz1);
tx1*=inv_d;
ty1*=inv_d;
tz1*=inv_d;
r=sqrt(r2);
rm1=(cnt+0)->radius-r;
if(r>1.e-6)
{
inv_d=1./r;
srd->x=((cnt+0)->radius*x1-rm1*(cnt+0)->cx)*inv_d;
srd->y=((cnt+0)->radius*y1-rm1*(cnt+0)->cy)*inv_d;
srd->z=((cnt+0)->radius*z1-rm1*(cnt+0)->cz)*inv_d;
}
else
{
r=1.e-6;
rp1=(cnt+0)->radius-r;
inv_d=1./r;
srd->x=((cnt+0)->radius*x1-rp1*(cnt+0)->cx)*inv_d;
srd->y=((cnt+0)->radius*y1-rp1*(cnt+0)->cy)*inv_d;
srd->z=((cnt+0)->radius*z1-rp1*(cnt+0)->cz)*inv_d;
}
velupdate(nx,ny,nz,tx1,ty1,tz1,a->n,srd->n);
return;
}
}
}
if(a->plane==length-1)
{
trvx=(cnt+a->plane)->dirx;
trvy=(cnt+a->plane)->diry;
trvz=(cnt+a->plane)->dirz;
dirvec=trvx*(a->cx-x1)+trvy*(a->cy-y1)+trvz*(a->cz-z1);
if(dirvec<0)
{
nx=(at+w2-1)->nx;
ny=(at+w2-1)->ny;
nz=(at+w2-1)->nz;
scalar=nx*((at+w2-1)->cx-x1)+ny*((at+w2-1)->cy-y1)+nz*((at+w2-1)->cz-z1);
srd->x=x1+nx*scalar;
srd->y=y1+ny*scalar;
srd->z=z1+nz*scalar;
velupdate2(a->n,srd->n);
return;
}
else
{
r2=dist2((cnt+length-1)->cx,x1,(cnt+length-1)->cy,y1,(cnt+length-1)->cz,z1);
if(r2<(cnt+length-1)->radius2)
{
nx=x1-(cnt+length-1)->cx;
ny=y1-(cnt+length-1)->cy;
nz=z1-(cnt+length-1)->cz;
inv_d=1./sqrt(nx*nx+ny*ny+nz*nz);
nx*=inv_d;
ny*=inv_d;
nz*=inv_d;
tx1=crossx(nx,ny,nz,((cnt+length-1)->cx-a->cx),((cnt+length-1)->cy-a->cy),((cnt+length-1)->cz-a->cz));
ty1=crossy(nx,ny,nz,((cnt+length-1)->cx-a->cx),((cnt+length-1)->cy-a->cy),((cnt+length-1)->cz-a->cz));
tz1=crossz(nx,ny,nz,((cnt+length-1)->cx-a->cx),((cnt+length-1)->cy-a->cy),((cnt+length-1)->cz-a->cz));
inv_d=1./sqrt(tx1*tx1+ty1*ty1+tz1*tz1);
tx1*=inv_d;
ty1*=inv_d;
tz1*=inv_d;
r=sqrt(r2);
rp1=(cnt+length-1)->radius-r;
if(r>1.e-6)
{
inv_d=1./r;
srd->x=((cnt+length-1)->radius*x1-rp1*(cnt+length-1)->cx)*inv_d;
srd->y=((cnt+length-1)->radius*y1-rp1*(cnt+length-1)->cy)*inv_d;
srd->z=((cnt+length-1)->radius*z1-rp1*(cnt+length-1)->cz)*inv_d;
}
else
{
r=1.e-6;
rp1=(cnt+length-1)->radius-r;
inv_d=1./r;
srd->x=((cnt+length-1)->radius*x1-rp1*(cnt+length-1)->cx)*inv_d;
srd->y=((cnt+length-1)->radius*y1-rp1*(cnt+length-1)->cy)*inv_d;
srd->z=((cnt+length-1)->radius*z1-rp1*(cnt+length-1)->cz)*inv_d;
}
velupdate(nx,ny,nz,tx1,ty1,tz1,a->n,srd->n);
return;
}
}
}
}
}
void restartb()
{
unsigned long int i,j,ic,ki,kj,pj;
double xl,yl,zl,eps;
struct redcell *rc;
static double *drhv; //////////////
static double *drhb; /////////////
static double *drhu; ////////////
static double *drha; ////////
static int n410st = 0;
if(n410st == 0){
drhv = (double *) malloc(sizeof(double)*(5000));//////////////
if(drhv == NULL) exit(1);////////////
drhb = (double *) malloc(sizeof(double)*(5000));///////////////
if(drhb == NULL) exit(1);////////////
drhu = (double *) malloc(sizeof(double)*(5000));///////////////
if(drhu == NULL) exit(1);/////////////
drha = (double *) malloc(sizeof(double)*(5000));///////////////
if(drha == NULL) exit(1);///////////
n410st++;
}
j=rank*500;
eps=0.0000000005;
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
rc->ip = 2;
}
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if((rc->xcm > xsl+Rrbc+0.5) && (rc->xcm < xer-Rrbc-0.5) )
{
rc->ip = 0;
rc->id=j;
j++;
}
}
kj=0;
drhv[kj++] = 0.0;
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->xcm >= xer-Rrbc-0.5)
{
rc->ip = 1;
rc->id=j;
drhv[kj++]=(double)rc->id;
drhv[kj++]=rc->xcm;
drhv[kj++]=rc->ycm;
drhv[kj++]=rc->zcm;
j++;
}
}
if(kj == 1) kj=0;
drhv[0]=(double)kj;
drhb[0] = 0.0;
if(numprocs > 1){
MPI_Sendrecv(&drhv[0],kj,MPI_DOUBLE,right,199,&drhb[0],5000,MPI_DOUBLE,left,199,comm1d,&status);
}
ki = (int)drhb[0];
kj = 1;
pj=0;
drhu[pj++] = 0.0;
while( kj < ki){
i=(int)drhb[kj++];
xl=drhb[kj++];
yl=drhb[kj++];
zl=drhb[kj++];
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if((rc->xcm > (xl-eps)) && (rc->xcm < (xl+eps)) && (rc->ycm > (yl-eps)) && (rc->ycm < (yl+eps)) && (rc->zcm > (zl-eps)) && (rc->zcm < (zl+eps)) && rc->xcm <= xsl+Rrbc+0.5)
{
rc->ip = 1;
rc->id = i;
ic=Nrbcp+10;
}
}
if(ic < Nrbcp+3)
{
drhu[pj++]=xl;
drhu[pj++]=yl;
drhu[pj++]=zl;
}
}
if(pj == 1) pj=0;
drhu[0]=(double)pj;
drha[0] = 0.0;
if(numprocs > 1){
MPI_Sendrecv(&drhu[0],pj,MPI_DOUBLE,left,201,&drha[0],5000,MPI_DOUBLE,right,201,comm1d,&status);
}
ki = (int)drha[0];
kj = 1;
while( kj < ki){
xl=drha[kj++];
yl=drha[kj++];
zl=drha[kj++];
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if((rc->xcm > (xl-eps)) && (rc->xcm < (xl+eps)) && (rc->ycm > (yl-eps)) && (rc->ycm < (yl+eps)) && (rc->zcm > (zl-eps)) && (rc->zcm < (zl+eps)))
{
rc->ip = 2;
}
}
}
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
rc->Vxcm=0.0;
rc->Vycm=0.0;
rc->Vzcm=0.0;
rc->wx=0.0;
rc->wy=0.0;
rc->wz=0.0;
rc->ex=1.0;
rc->ey=0.0;
rc->ez=0.0;
rc->Momx=0.0;
rc->Momy=0.0;
rc->Momz=0.0;
rc->Amomx=0.0;
rc->Amomy=0.0;
rc->Amomz=0.0;
printf(" %i\t%i\t%i\t%i\t%i\t%le\t%le\t%le\n",rank,Nrbcp,ic,rc->id,rc->ip,rc->xcm,rc->ycm,rc->zcm);
}
}
void restart()
{
unsigned long int i,j,ic,kj,typstart,ip,id;
double xstart,ystart,zstart,vxstart,vystart,vzstart,axstart,aystart,azstart,xcmstart,ycmstart,zcmstart,Vxcmstart,Vycmstart,Vzcmstart,wxstart,wystart,wzstart,exstart,eystart,ezstart;
struct aftr *a,*a1;
struct sol *srd;
struct redcell *rc;
FILE *fp;
char readfile[256];
sprintf(readfile,"%s/sol_atp%ld-%09ld.dat",fname,rank,tsim);
printf("rank=%d readfile is=%s\n",rank,readfile);
fp = fopen(readfile,"r");
if(fp == NULL){
printf("Read file does not exist.");
exit(1);
}
fscanf(fp,"%le",&(totdt));
for(i=0;i<Naf;i++)
{
fscanf(fp,"%ld\t%ld\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le",&(j),&(typstart),&(xstart),&(ystart),&(zstart),&(vxstart),&(vystart),&(vzstart),&(axstart),&(aystart),&(azstart));
if(i != j) printf("RestartERRRRRRRRRRRRRRRRRRRROOOOOOOOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRR");
a=at+i;
a->typ=typstart;
a->xreal=xstart;
a->yreal=ystart;
a->zreal=zstart;
a->vx=vxstart;
a->vy=vystart;
a->vz=vzstart;
a->ax=axstart;
a->ay=aystart;
a->az=azstart;
}
fscanf(fp,"%ld",&(Nrbcp));
for(ic=0;ic<Nrbcp;ic++)
{
// fscanf(fp,"%ld\t%le\t%le\t%le",&(j),&(xcmstart),&(ycmstart),&(zcmstart));
fscanf(fp,"%ld\t%ld\t%ld\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le",&(j),&(ip),&(id),&(xcmstart),&(ycmstart),&(zcmstart),&(Vxcmstart),&(Vycmstart),&(Vzcmstart),&(wxstart),&(wystart),&(wzstart),&(exstart),&(eystart),&(ezstart));
if(ic != j) printf("RestartERRRRRRRRRRRRRRRRRRRROOOOOOOOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRR1");
rc=rbc+ic;
rc->xcm=xcmstart;
rc->ycm=ycmstart;
rc->zcm=zcmstart;
rc->id=id;
rc->ip=ip;
rc->Vxcm=Vxcmstart;
rc->Vycm=Vycmstart;
rc->Vzcm=Vzcmstart;
rc->wx=wxstart;
rc->wy=wystart;
rc->wz=wzstart;
rc->ex=exstart;
rc->ey=eystart;
rc->ez=ezstart;
rc->Momx=0.0;
rc->Momy=0.0;
rc->Momz=0.0;
rc->Amomx=0.0;
rc->Amomy=0.0;
rc->Amomz=0.0;
}
fscanf(fp,"%ld",&(Nsolp));
for(i=0;i<Nsolp;i++)
{
fscanf(fp,"%ld\t%ld\t%le\t%le\t%le\t%le\t%le\t%le",&(j),&(kj),&(xstart),&(ystart),&(zstart),&(vxstart),&(vystart),&(vzstart));
if(i != j || i != kj) printf("RestartERRRRRRRRRRRRRRRRRRRROOOOOOOOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRR2");
srd=ss+i;
srd->n=kj;
srd->x=xstart;
srd->y=ystart;
srd->z=zstart;
srd->vx=vxstart;
srd->vy=vystart;
srd->vz=vzstart;
}
fclose(fp);
}
void snapshot1()
{
unsigned long int i,ic,j,k,j1,j2,m,out,nafp=0;
double xl,yl,zl,norm,x1,y1,z1,x2,y2,z2,x3,y3,z3,inv_d,d,norm1,d2;
double tx,ty,tz;
struct aftr *a,*a1;
struct cnline *c;
struct sol *srd;
struct redcell *rc;
FILE *fp;
char outfile[256];
j=length-1;
sprintf(outfile,"%s/sol_atp%ld-%09ld.dat",foutpoint1,rank,tsim);
fp = fopen(outfile,"w");
if(fp == NULL)
{
printf("EEEEEEEERRRRRRRRRRRRROOOOOOOOORRRRRRRR");
exit(1);
}
fprintf(fp,"%le\n",totdt);
for(i=0;i<Naf;i++)
{
a=at+i;
fprintf(fp,"%ld\t%ld\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\n",i,a->typ,a->xreal,a->yreal,a->zreal,a->vx,a->vy,a->vz,a->ax,a->ay,a->az);
}
fprintf(fp,"%ld\n",Nrbcp);
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
fprintf(fp,"%ld\t%ld\t%ld\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\n",ic,rc->ip,rc->id,rc->xcm,rc->ycm,rc->zcm,rc->Vxcm,rc->Vycm,rc->Vzcm,rc->wx,rc->wy,rc->wz,rc->ex,rc->ey,rc->ez);
}
fprintf(fp,"%ld\n",Nsolp);
for(i=0;i<Nsolp;i++)
{
srd=ss+i;
fprintf(fp,"%ld\t%ld\t%le\t%le\t%le\t%le\t%le\t%le\n",i,srd->n,srd->x,srd->y,srd->z,srd->vx,srd->vy,srd->vz);
}
fclose(fp);
}
void snapshot()
{
unsigned long int i,ic,j,k,j1,j2,m,out,nafp=0;
double xl,yl,zl,norm,x1,y1,z1,x2,y2,z2,x3,y3,z3,inv_d,d,norm1,d2;
double tx,ty,tz;
struct aftr *a,*a1;
struct cnline *c;
struct sol *srd;
struct redcell *rc;
FILE *fp;
char outfile[256];
j=length-1;
sprintf(outfile,"%s/atpic%ld-%09ld.dat",foutpoint,rank,tsim);
fp = fopen(outfile,"w");
if(fp == NULL)
{
printf("EEEEEEEERRRRRRRRRRRRROOOOOOOOORRRRRRRR");
exit(1);
}
for(i=0;i<Naf;i++)
{
a=at+i;
if(a->typ == 0) nafp++;
}
fprintf(fp,"%ld\n",nafp);
for(i=0;i<Naf-nflagella;i++)
{
a=at+i;
if(a->typ != 0) continue;
fprintf(fp,"%ld\t%le\t%le\t%le\n",i,a->xreal,a->yreal,a->zreal);
/* if(i==j)
{
fprintf(fp,"\n");
j+=length;
}
*/
}
m=0;
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
if((rc->xcm <= xer) && (rc->xcm > xsl))
{
m++;
}
}
fprintf(fp,"%ld\n",m);
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
if((rc->xcm <= xer) && (rc->xcm > xsl))
{
fprintf(fp,"%le\t%le\t%le\n",rc->xcm,rc->ycm,rc->zcm);
}
}
fclose(fp);
/*
for(i=0;i<length;i++)
{
fprintf(fp,"%le\t%le\t%le\n",(at+i)->xreal,(at+i)->yreal,(at+i)->zreal);
}
fprintf(fp,"\n\n");
a=(at+1);
i=1;
do
{
if(i==1)
a1=a;
if(i==2)
a1=a->helix[0];
if(i>2)
a1=a->helix[1];
fprintf(fp,"%le\t%le\t%le\n",a1->xreal,a1->yreal,a1->zreal);
//printf("%ld\t%ld\n",a->n,a1->n);
a=a1;
i++;
}
while(i<length+nflagella-5);
fprintf(fp,"\n\n");
for(i=0;i<flalen;i++)
{
c=cnt+i;
fprintf(fp,"%le\t%le\t%le\n",c->xreal,c->yreal,c->zreal);
} */
/*
j=length-1;
sprintf(name,"c-%09ld.dat",tsim);
fp=fopen(name,"w");
for(i=0;i<Naf-nflagella;i++)
{
a=at+i;
fprintf(fp,"%le\t%le\t%le\n",a->cx,a->cy,a->cz);
if(i==j)
{
fprintf(fp,"\n");
j+=length;
}
}
for(i=0;i<length;i++)
{
fprintf(fp,"%le\t%le\t%le\n",(at+i)->cx,(at+i)->cy,(at+i)->cz);
}
fprintf(fp,"\n\n");
a=(at+1);
i=1;
do
{
if(i==1)
a1=a;
if(i==2)
a1=a->helix[0];
if(i>2)
a1=a->helix[1];
fprintf(fp,"%le\t%le\t%le\n",a1->cx,a1->cy,a1->cz);
//printf("%ld\t%ld\n",a->n,a1->n);
a=a1;
i++;
}
while(i<length+nflagella-5);
fclose(fp);
*/
/*j=length-1;
sprintf(name,"c-%09ld.dat",tsim);
fp=fopen(name,"w");
for(i=0;i<Naf;i++)
{
a=at+i;
fprintf(fp,"%le\t%le\t%le\n",a->cx,a->cy,a->cz);
if(i==j)
{
fprintf(fp,"\n");
j+=length;
}
}
for(i=0;i<length;i++)
{
fprintf(fp,"%le\t%le\t%le\n",(at+i)->cx,(at+i)->cy,(at+i)->cz);
}
fprintf(fp,"\n\n");
for(i=0;i<flalen;i++)
{
c=cnt+i;
fprintf(fp,"%le\t%le\t%le\n",c->xreal,c->yreal,c->zreal);
}
fclose(fp);
*/
/* sprintf(name,"velacc-%09ld.dat",tsim);
fp=fopen(name,"w");
a=(at+1);
i=1;
do
{
if(i==1)
a1=a;
if(i==2)
a1=a->helix[0];
if(i>2)
a1=a->helix[1];
fprintf(fp,"%ld\t%le\t%le\t%le\t%le\t%le\t%le\n",i,a1->vx,a1->vy,a1->vz,a1->ax,a1->ay,a1->az);
//printf("%ld\t%ld\n",a->n,a1->n);
a=a1;
i++;
}
while(i<length+nflagella-5);
fprintf(fp,"\n\n");
tx=ty=tz=0;
for(i=0;i<Naf;i++)
{
a=at+i;
fprintf(fp,"%ld\t%le\t%le\t%le\n",i,a->vx,a->vy,a->vz);
}
fprintf(fp,"\n\n");
for(i=0;i<Naf;i++)
{
a=at+i;
fprintf(fp,"%ld\t%le\t%le\t%le\n",i,a->ax,a->ay,a->az);
}
fprintf(fp,"\n\n");
for(i=0;i<Naf;i++)
{
a=at+i;
fprintf(fp,"%ld\t%le\t%le\t%le\n",i,crossx(a->xreal,a->yreal,a->zreal,a->ax,a->ay,a->az),crossy(a->xreal,a->yreal,a->zreal,a->ax,a->ay,a->az),crossz(a->xreal,a->yreal,a->zreal,a->ax,a->ay,a->az));
}
fprintf(fp,"\n\n");
for(i=0;i<Naf;i++)
{
a=at+i;
inv_d=1./(a->xreal*a->xreal+a->yreal*a->yreal+a->zreal*a->zreal);
fprintf(fp,"%ld\t%le\t%le\t%le\n",i,crossx(a->xreal,a->yreal,a->zreal,a->vx,a->vy,a->vz)*inv_d,crossy(a->xreal,a->yreal,a->zreal,a->vx,a->vy,a->vz)*inv_d,crossz(a->xreal,a->yreal,a->zreal,a->vx,a->vy,a->vz)*inv_d);
}
fclose(fp);
*/
/*
sprintf(name,"sol-%09ld.dat",tsim);
fp=fopen(name,"w");
for(k=0;k<Nsol;k++)
{
srd=ss+k;
//srd->sp=0;
for(i=0;i<length;i++)
{
c=cnt+i;
//for(j=0;j<8;j++)
//{
//d2=dist2(c->cx,srd->pbcx[j],c->cy,srd->pbcy[j],c->cz,srd->pbcz[j]);
d2=dist2(c->cx,srd->x,c->cy,srd->y,c->cz,srd->z);
if(d2<5.)
{
x1=srd->x;
y1=srd->y;
z1=srd->z;
fprintf(fp,"%le\t%le\t%le\n",x1,y1,z1);
//srd->sp=1;
//}
}
}
}
fclose(fp);
*/
/*
if(tsim%20000==0 && tsim>0)
{
sprintf(name,"flow-%09ld.dat",tsim);
fp=fopen(name,"w");
for(i=0;i<L;i++)
for(j=0;j<L;j++)
for(k=0;k<L;k++)
{
tx=N_VCM[i+L*(j+L*k)].Vx;
ty=N_VCM[i+L*(j+L*k)].Vy;
tz=N_VCM[i+L*(j+L*k)].Vz;
d=sqrt(tx*tx+ty*ty+tz*tz);
if(d>0)
{
inv_d=1./d;
tx*=inv_d;
ty*=inv_d;
tz*=inv_d;
fprintf(fp,"%le\t%le\t%le\t%le\t%le\t%le\n",(double)i+0.5,(double)j+0.5,(double)k+0.5,tx,ty,tz);
//printf("%le\t%le\t%le\t%le\t%le\t%le\n",(double)i,(double)L*(double)i,(double)L*(double)L*(double)i,tx,ty,tz);
}
}
fclose(fp);
//exit(1);
}*/
}
void radiuscal()
{
unsigned long int i,j,k;
double dam1x,dam1y,dam1z,fx,fy,fz,l0,d2,d,inv_d,dax,day,daz,df,ds;
struct aftr *a1,*a2,*a3,*a4;
for(i=0;i<length-1;i++)
{
k=i;
j=0;
do
{
a1=(at+(k+j*length));
a2=(at+k+length*(j+1)+1);
a3=(at+(k+j*length+length));
a4=(at+(k+length*(j+1)+1)-length);
l0=a1->bn;
dam1x=(a1->xreal-a2->xreal);
dam1y=(a1->yreal-a2->yreal);
dam1z=(a1->zreal-a2->zreal);
d2=dist2(a1->xreal,a2->xreal,a1->yreal,a2->yreal,a1->zreal,a2->zreal);
df=sqrt(d2);
a1->rad2=d2*0.25;
a1->rad=df*0.5;
a2->rad2=d2*0.25;
a2->rad=df*0.5;
dax=(a3->xreal-a4->xreal);
day=(a3->yreal-a4->yreal);
daz=(a3->zreal-a4->zreal);
d2=dist2(a3->xreal,a4->xreal,a3->yreal,a4->yreal,a3->zreal,a4->zreal);
ds=sqrt(d2);
a3->rad2=d2*0.25;
a3->rad=ds*0.5;
a4->rad2=d2*0.25;
a4->rad=ds*0.5;
j++;
//printf("%ld\t%ld\t%ld\t%ld\n",a1->n,a2->n,a3->n,a4->n);
}
while(j<npoints-1);
}
//exit(1);
}
void bendbodyforcet0()
{
unsigned long int i,j,k;
double mrix,mriy,mriz,mrip1x,mrip1y,mrip1z;
double rix,riy,riz,rip1x,rip1y,rip1z,inv_ri,inv_rip1,inv_mri,inv_mrip1;
double inv_delta,cp,cm,delta;
double fx,fy,fz,tfx,tfy,tfz,ttx,tty,ttz,Tx,Ty,Tz;
struct aftr *a1,*a2,*a3;
j=0;
do
{
for(i=0;i<length-2;i++)
{
a1=at+j*length+i;
a2=at+j*length+i+1;
a3=at+j*length+i+2;
rix=a2->xreal-a1->xreal;
riy=a2->yreal-a1->yreal;
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal-a2->zreal;
/*
rix=a2->x0-a1->x0;
riy=a2->y0-a1->y0;
riz=a2->z0-a1->z0;
rip1x=a3->x0-a2->x0;
rip1y=a3->y0-a2->y0;
rip1z=a3->z0-a2->z0;
*/
inv_ri=1./sqrt(rix*rix+riy*riy+riz*riz);
inv_rip1=1./sqrt(rip1x*rip1x+rip1y*rip1y+rip1z*rip1z);
rix*=inv_ri;
riy*=inv_ri;
riz*=inv_ri;
rip1x*=inv_rip1;
rip1y*=inv_rip1;
rip1z*=inv_rip1;
cp=rix*rip1x+riy*rip1y+riz*rip1z;
a2->cangle1=cp;
// if(rank==1) printf("%ld\t%le\n",a2->n,acos(cp)*180./PI);
}
//printf("\n\n");
j++;
}
while(j<npoints);
//exit(1);
}
void ND()
{
unsigned long int i,ic,j,ki=0,kj=0,pi=0,pj=0,tsi=0,tsj=0,icol,jcol,iout;
double xl,yl,zl,sx,sy,sz,nfx,nfy,nfz,rbp;
double xcm,ycm,zcm,vxcm,vycm,vzcm,fx,fy,fz;
double xncm,yncm,zncm,rpc;
double xxcm,yycm,zzcm,vrx,vry,vrz,bici,deltaa,dtnew,rsq,vsq,tcol,tcol2,tcolmin,vradi,vdir,dryp,dic,dtnu;
struct aftr *a;
struct redcell *rc;
struct redcell *rc2;
nfx=nfy=nfz=1./(double)tsim;
xcm=ycm=zcm=vxcm=vycm=vzcm=xncm=yncm=zncm=0;
fx=fy=fz=j=0;
//////////////////////////
static double *hu;
static double *hv;
static double *ha;
static double *hb;
static double *pshu;
static double *pshv;
static double *psha;
static double *pshb;
static double *tshu;
static double *tshv;
static double *tsha;
static double *tshb;
static int n1st = 0;
if(n1st == 0){
hu = (double *) malloc(sizeof(double)*(5000));
if(hu == NULL){
printf("\n err_message(1)");
exit(1);
}
hv = (double *) malloc(sizeof(double)*(5000));
if(hv == NULL) exit(1);
ha = (double *) malloc(sizeof(double)*(5000));
if(ha == NULL) exit(1);
hb = (double *) malloc(sizeof(double)*(5000));
if(hb == NULL) exit(1);
pshu = (double *) malloc(sizeof(double)*(5000));
if(pshu == NULL){
printf("\n err_message(1)");
exit(1);
}
pshv = (double *) malloc(sizeof(double)*(5000));
if(pshv == NULL) exit(1);
psha = (double *) malloc(sizeof(double)*(5000));
if(psha == NULL) exit(1);
pshb = (double *) malloc(sizeof(double)*(5000));
if(pshb == NULL) exit(1);
tshu = (double *) malloc(sizeof(double)*(5000));
if(tshu == NULL){
printf("\n err_message(1)");
exit(1);
}
tshv = (double *) malloc(sizeof(double)*(5000));
if(tshv == NULL) exit(1);
tsha = (double *) malloc(sizeof(double)*(5000));
if(tsha == NULL) exit(1);
tshb = (double *) malloc(sizeof(double)*(5000));
if(tshb == NULL) exit(1);
n1st++;
}
// MPI_Status status;
pshu[pi++] = 0.0;
pshv[pj++] = 0.0;
psha[0] = 0.0;
pshb[0] = 0.0;
tshu[tsi++] = 0.0;
tshv[tsj++] = 0.0;
tsha[0] = 0.0;
tshb[0] = 0.0;
hu[ki++] = 0.0;
hv[kj++] = 0.0;
ha[0] = 0.0;
hb[0] = 0.0;
/////////////////////////////////////////
for(i=0;i<Naf;i++)
{
a=at+i;
///////////////////////////////////
if(a->typ != 0)continue;
a->vx+=(0.5*dtaf*a->ax);
a->vy+=(0.5*dtaf*a->ay);
a->vz+=(0.5*dtaf*a->az);
sx=dtaf*a->vx;
sy=dtaf*a->vy;
sz=dtaf*a->vz;
a->xreal+=sx;
a->yreal+=sy;
a->zreal+=sz;
xl=a->xreal;
yl=a->yreal;
zl=a->zreal;
rbp=sqrt((yl-ylc)*(yl-ylc)+(zl-zlc)*(zl-zlc));
if(rbp > Lr0)
{
bounceback(1,i,dtaf);
}
xl=a->xreal;
yl=a->yreal;
zl=a->zreal;
if(xl>(double)L || xl<0 || yl>(double)Lr1 || yl<0 || zl>(double)Lr1 || zl<0)
{
fprintf(fp,"cellbody out of the box\t%le\t%le\t%le\n",xl,yl,zl);
fp=fopen("info.txt","a");
fprintf(fp,"cellbody out of the box\t%le\t%le\t%le\n",xl,yl,zl);
fclose(fp);
exit(1);
}
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
rpc=sqrt((xl-rc->xcm)*(xl-rc->xcm)+(yl-rc->ycm)*(yl-rc->ycm)+(zl-rc->zcm)*(zl-rc->zcm));
if(rpc<Rrbc)
{
// printf(" rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrraaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaannnnnnnnnnnnnnnnnnnnnnnkkkkkkkkkkkkkkkkkkkk=%i\n",rank);
bouncebackrbc(1,i,ic,dtaf);
ic=50000;
}
}
////////////////////
xl=a->xreal;
yl=a->yreal;
zl=a->zreal;
///////////////////////////////
tshv[tsj++]=(double)i;
tshv[tsj++]=a->xreal;
tshv[tsj++]=a->yreal;
tshv[tsj++]=a->zreal;
tshu[tsi++]=(double)i;
tshu[tsi++]=a->xreal;
tshu[tsi++]=a->yreal;
tshu[tsi++]=a->zreal;
if(xl >= xer)
{
a->typ=1;
pshv[pj++]=(double)i;
pshv[pj++]=a->xreal;
pshv[pj++]=a->yreal;
pshv[pj++]=a->zreal;
pshv[pj++]=a->vx;
pshv[pj++]=a->vy;
pshv[pj++]=a->vz;
}
else
{
hv[kj++]=(double)i;
hv[kj++]=a->xreal;
hv[kj++]=a->yreal;
hv[kj++]=a->zreal;
}
if(xl < xsl)
{
a->typ=1;
pshu[pi++]=(double)i;
pshu[pi++]=a->xreal;
pshu[pi++]=a->yreal;
pshu[pi++]=a->zreal;
pshu[pi++]=a->vx;
pshu[pi++]=a->vy;
pshu[pi++]=a->vz;
}
else
{
hu[ki++]=(double)i;
hu[ki++]=a->xreal;
hu[ki++]=a->yreal;
hu[ki++]=a->zreal;
}
}
if(ki == 1) ki=0;
if(kj == 1) kj=0;
hu[0]=(double)ki;
hv[0]=(double)kj;
if(pi == 1) pi=0;
if(pj == 1) pj=0;
pshu[0]=(double)pi;
pshv[0]=(double)pj;
if(numprocs > 1){
MPI_Sendrecv(&hu[0],ki,MPI_DOUBLE,left,1,&ha[0],5000,MPI_DOUBLE,right,1,comm1d,&status);
MPI_Sendrecv(&hv[0],kj,MPI_DOUBLE,right,2,&hb[0],5000,MPI_DOUBLE,left,2,comm1d,&status);
MPI_Sendrecv(&pshu[0],pi,MPI_DOUBLE,left,100,&psha[0],5000,MPI_DOUBLE,right,100,comm1d,&status);
MPI_Sendrecv(&pshv[0],pj,MPI_DOUBLE,right,101,&pshb[0],5000,MPI_DOUBLE,left,101,comm1d,&status);
}
///*****************************************//
ki = (int)ha[0];
kj = 1;
while( kj < ki){
i=(int)ha[kj++];
a=at+i;
a->typ=1;
a->xreal=ha[kj++];
a->yreal=ha[kj++];
a->zreal=ha[kj++];
// printf(" %i\t%i\t%i\t%le\t%i\t%i\n",ki,rank,a->typ,a->xreal,xsl,xer);
}
ki = (int)hb[0];
kj = 1;
while( kj < ki){
i=(int)hb[kj++];
a=at+i;
a->typ=1;
a->xreal=hb[kj++];
a->yreal=hb[kj++];
a->zreal=hb[kj++];
}
///***************************************//
ki = (int)psha[0];
kj = 1;
while( kj < ki){
i=(int)psha[kj++];
a=at+i;
a->typ=0;
a->xreal=psha[kj++];
a->yreal=psha[kj++];
a->zreal=psha[kj++];
a->vx=psha[kj++];
a->vy=psha[kj++];
a->vz=psha[kj++];
// printf(" %i\t%i\t%i\t%le\t%i\t%i\n",ki,rank,a->typ,a->xreal,xsl,xer);
}
ki = (int)pshb[0];
kj = 1;
while( kj < ki){
i=(int)pshb[kj++];
a=at+i;
a->typ=0;
a->xreal=pshb[kj++];
a->yreal=pshb[kj++];
a->zreal=pshb[kj++];
a->vx=pshb[kj++];
a->vy=pshb[kj++];
a->vz=pshb[kj++];
// printf(" %i\t%i\t%i\t%le\t%i\t%i\n",ki,rank,a->typ,a->xreal,xsl,xer);
}
///*****************************************///
if(tsi == 1) tsi=0;
if(tsj == 1) tsj=0;
tshu[0]=(double)tsi;
tshv[0]=(double)tsj;
if(rank > 1) MPI_Send(&tshu[0],tsi,MPI_DOUBLE,rank-2,102,MPI_COMM_WORLD);
if(rank < numprocs-2)
{
MPI_Recv(&tsha[0],5000,MPI_DOUBLE,rank+2,102,MPI_COMM_WORLD,&status);
ki = (int)tsha[0];
kj = 1;
while( kj < ki){
i=(int)tsha[kj++];
a=at+i;
a->typ=1;
a->xreal=tsha[kj++];
a->yreal=tsha[kj++];
a->zreal=tsha[kj++];
// printf(" %i\t%i\t%i\t%le\t%i\t%i\n",ki,rank,a->typ,a->xreal,xsl,xer);
}
}
if(rank < numprocs-2) MPI_Send(&tshv[0],tsj,MPI_DOUBLE,rank+2,103,MPI_COMM_WORLD);
if(rank > 1)
{
MPI_Recv(&tshb[0],5000,MPI_DOUBLE,rank-2,103,MPI_COMM_WORLD,&status);
ki = (int)tshb[0];
kj = 1;
while( kj < ki){
i=(int)tshb[kj++];
a=at+i;
a->typ=1;
a->xreal=tshb[kj++];
a->yreal=tshb[kj++];
a->zreal=tshb[kj++];
// printf(" %i\t%i\t%i\t%le\t%i\t%i\n",ki,rank,a->typ,a->xreal,xsl,xer);
}
}
//***********************************************//
calculateforces();
}
void MDRBC()
{
unsigned long int i,ic,j,ki=0,kj=0,pi=0,pj=0,tsi=0,tsj=0,icol,jcol,iout,out;
double xl,yl,zl,sx,sy,sz,nfx,nfy,nfz,rbp,rbpo,xlo,ylo,zlo,dtp;
double xcm,ycm,zcm,vxcm,vycm,vzcm,fx,fy,fz;
double xncm,yncm,zncm,rpc,rctoc;
double xxcm,yycm,zzcm,vrx,vry,vrz,bici,deltaa,dtnew,rsq,vsq,tcol,tcol2,tcolmin,vradi,vdir,dryp,dic,dtnu;
struct aftr *a;
struct redcell *rc;
struct redcell *rc2;
nfx=nfy=nfz=1./(double)tsim;
xcm=ycm=zcm=vxcm=vycm=vzcm=xncm=yncm=zncm=0;
fx=fy=fz=j=0;
/////////////////////////////
static double *rhu;
static double *rhv;
static double *rha;
static double *rhb;
static double *rpshu;
static double *rpshv;
static double *rpsha;
static double *rpshb;
static int n310st = 0;
if(n310st == 0){
rhu = (double *) malloc(sizeof(double)*(5000));
if(rhu == NULL){
printf("\n err_message(1)");
exit(1);
}
rhv = (double *) malloc(sizeof(double)*(5000));
if(rhv == NULL) exit(1);
rha = (double *) malloc(sizeof(double)*(5000));
if(rha == NULL) exit(1);
rhb = (double *) malloc(sizeof(double)*(5000));
if(rhb == NULL) exit(1);
rpshu = (double *) malloc(sizeof(double)*(5000));
if(rpshu == NULL){
printf("\n err_message(1)");
exit(1);
}
rpshv = (double *) malloc(sizeof(double)*(5000));
if(rpshv == NULL) exit(1);
rpsha = (double *) malloc(sizeof(double)*(5000));
if(rpsha == NULL) exit(1);
rpshb = (double *) malloc(sizeof(double)*(5000));
if(rpshb == NULL) exit(1);
n310st++;
}
//**********************************************//
//////////////////////////////
if(OBSM==1)
{
ki=0;
kj=0;
rhu[ki++] = 0.0;
rhv[kj++] = 0.0;
rha[0] = 0.0;
rhb[0] = 0.0;
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip != 1)continue;
if(rc->xcm <= xsl+Rrbc+0.55)
{
rhu[ki++]=(double)rc->id;
rhu[ki++]=rc->Momx;
rhu[ki++]=rc->Momy;
rhu[ki++]=rc->Momz;
rhu[ki++]=rc->Amomx;
rhu[ki++]=rc->Amomy;
rhu[ki++]=rc->Amomz;
rhu[ki++]=rc->Vxcm;
rhu[ki++]=rc->Vycm;
rhu[ki++]=rc->Vzcm;
rhu[ki++]=rc->wx;
rhu[ki++]=rc->wy;
rhu[ki++]=rc->wz;
}
else if(rc->xcm >= xer-Rrbc-0.55)
{
rhv[kj++]=(double)rc->id;
rhv[kj++]=rc->Momx;
rhv[kj++]=rc->Momy;
rhv[kj++]=rc->Momz;
rhv[kj++]=rc->Amomx;
rhv[kj++]=rc->Amomy;
rhv[kj++]=rc->Amomz;
}
}
if(ki == 1) ki=0;
if(kj == 1) kj=0;
rhu[0]=(double)ki;
rhv[0]=(double)kj;
if(numprocs > 1){
MPI_Sendrecv(&rhu[0],ki,MPI_DOUBLE,left,185,&rha[0],5000,MPI_DOUBLE,right,185,comm1d,&status);
MPI_Sendrecv(&rhv[0],kj,MPI_DOUBLE,right,186,&rhb[0],5000,MPI_DOUBLE,left,186,comm1d,&status);
}
///*****************************************//
ki = (int)rha[0];
kj = 1;
while( kj < ki){
i=(int)rha[kj++];
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip != 1)continue;
if(rc->id == i)
{
rc->Momx+=rha[kj++];
rc->Momy+=rha[kj++];
rc->Momz+=rha[kj++];
rc->Amomx+=rha[kj++];
rc->Amomy+=rha[kj++];
rc->Amomz+=rha[kj++];
rc->Vxcm=rha[kj++];
rc->Vycm=rha[kj++];
rc->Vzcm=rha[kj++];
rc->wx=rha[kj++];
rc->wy=rha[kj++];
rc->wz=rha[kj++];
ic=Nrbcp+10;
}
}
if(ic < Nrbcp+3)
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_ip");
printf(" rank=%i id=%i xsl=%i xer=%i \n",rank,i,xsl,xer);
exit(1);
}
}
ki = (int)rhb[0];
kj = 1;
while( kj < ki){
i=(int)rhb[kj++];
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip != 1)continue;
if(rc->id == i)
{
rc->Momx+=rhb[kj++];
rc->Momy+=rhb[kj++];
rc->Momz+=rhb[kj++];
rc->Amomx+=rhb[kj++];
rc->Amomy+=rhb[kj++];
rc->Amomz+=rhb[kj++];
ic=Nrbcp+10;
}
}
if(ic < Nrbcp+3)
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_ipbb");
printf(" rank=%i id=%i xsl=%i xer=%i \n",rank,i,xsl,xer);
exit(1);
}
}
///***************************************//
/////////////////////
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
rc->Vxcm+=-(rc->Momx/Mrbc);
rc->Vycm+=-(rc->Momy/Mrbc);
rc->Vzcm+=-(rc->Momz/Mrbc);
rc->wx+=-(rc->Amomx/Irbc);
rc->wy+=-(rc->Amomy/Irbc);
rc->wz+=-(rc->Amomz/Irbc);
if(rc->Vxcm >1.0)rc->Vxcm=1.0;
if(rc->Vycm >1.0)rc->Vycm=1.0;
if(rc->Vzcm >1.0)rc->Vzcm=1.0;
if(rc->wx >1.0)rc->wx=1.0;
if(rc->wy >1.0)rc->wy=1.0;
if(rc->wz >1.0)rc->wz=1.0;
// rc->xcm+=rc->Vxcm*dt;
// rc->ycm+=rc->Vycm*dt;
// rc->zcm+=rc->Vzcm*dt;
rc->ex+=(crossx(rc->wx,rc->wy,rc->wz,rc->ex,rc->ey,rc->ez))*dt;
rc->ey+=(crossy(rc->wx,rc->wy,rc->wz,rc->ex,rc->ey,rc->ez))*dt;
rc->ez+=(crossz(rc->wx,rc->wy,rc->wz,rc->ex,rc->ey,rc->ez))*dt;
rc->ex=rc->ex/sqrt(rc->ex*rc->ex+rc->ey*rc->ey+rc->ez*rc->ez);
rc->ey=rc->ey/sqrt(rc->ex*rc->ex+rc->ey*rc->ey+rc->ez*rc->ez);
rc->ez=rc->ez/sqrt(rc->ex*rc->ex+rc->ey*rc->ey+rc->ez*rc->ez);
rc->Momx=0.0;
rc->Momy=0.0;
rc->Momz=0.0;
rc->Amomx=0.0;
rc->Amomy=0.0;
rc->Amomz=0.0;
}
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
rc->coltim=100.0;
}
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
rbp=sqrt((ylc-rc->ycm)*(ylc-rc->ycm)+(zlc-rc->zcm)*(zlc-rc->zcm));
if(rbp+Rrbc >= Lr0)
{
printf("aEsRRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_rbp");
printf(" rank=%i ic=%i ip=%i rbp=%e \n",rank,ic,rc->ip,rbp);
exit(1);
}
rbp=sqrt((ylc-(rc->ycm+rc->Vycm*dt))*(ylc-(rc->ycm+rc->Vycm*dt))+(zlc-(rc->zcm+rc->Vzcm*dt))*(zlc-(rc->zcm+rc->Vzcm*dt)));
if(rbp+Rrbc >= Lr0)
{
yycm=rc->ycm-ylc;
zzcm=rc->zcm-zlc;
vry=rc->Vycm;
vrz=rc->Vzcm;
bici=yycm*vry+zzcm*vrz;
if(bici >= 0.0)
{
rsq=yycm*yycm+zzcm*zzcm;
vsq=vry*vry+vrz*vrz;
dic=Lr0-Rrbc;
deltaa=bici*bici-vsq*(rsq-dic*dic);
if(deltaa>0)
{
tcol=(-bici+sqrt(deltaa))/vsq;
if(tcol < rc->coltim && tcol>0)
{
rc->coltim=tcol;
rc->partnr=5000000;
}
}
if(deltaa<=0.0 || tcol<=0.0)
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_deltat or tcol");
printf(" rank=%i ic=%i xsl=%i xer=%i deltaa=%e tcol=%e \n",rank,ic,xsl,xer,deltaa,tcol);
exit(1);
}
}
else
{
printf("ddddERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_deltat or tcol");
printf(" rank=%i ic=%i bici=%e tsim=%d \n",rank,ic,bici,tsim);
exit(1);
}
}
}
////////////////////////////////
for(ic=0;ic<Nrbcp-1;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
for(i=ic+1;i<Nrbcp;i++)
{
rc2=rbc+i;
if(rc2->ip == 2)continue;
xxcm=rc->xcm-rc2->xcm;
yycm=rc->ycm-rc2->ycm;
zzcm=rc->zcm-rc2->zcm;
vrx=rc->Vxcm-rc2->Vxcm;
vry=rc->Vycm-rc2->Vycm;
vrz=rc->Vzcm-rc2->Vzcm;
bici=xxcm*vrx+yycm*vry+zzcm*vrz;
if(bici < 0.0)
{
rsq=xxcm*xxcm+yycm*yycm+zzcm*zzcm;
vsq=vrx*vrx+vry*vry+vrz*vrz;
deltaa=bici*bici-vsq*(rsq-4*Rrbc*Rrbc);
if(deltaa>0)
{
tcol=(-bici-sqrt(deltaa))/vsq;
if(tcol < rc->coltim && tcol>0)
{
rc->coltim=tcol;
rc->partnr=i;
}
if(tcol < rc2->coltim && tcol>0)
{
rc2->coltim=tcol;
rc2->partnr=ic;
}
}
}
}
}
//////////////////////////
tcol=100.0;
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
if(rc->coltim < tcol)
{
tcol=rc->coltim;
icol=ic;
}
}
// rc=rbc+icol;
// jcol=rc->partnr;
// MPI_Allreduce(&a,&a,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
MPI_Allreduce(&tcol,&tcolmin,1,MPI_DOUBLE,MPI_MIN,comm1d);
//////////////////////////////////////////////////////////////////////
dtnew=dt;
dtnu=dt;
iout=1;
if(tcolmin <dt)
{
dtnew=tcolmin;
iout=0;
}
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
rc->xcm+=rc->Vxcm*dtnew;
rc->ycm+=rc->Vycm*dtnew;
rc->zcm+=rc->Vzcm*dtnew;
rc->coltim-=dtnew;
}
// if(rc->id == 369 && tsim >459000)
//////////////////////////////////////////////////////////////////////////////////////////
while(iout==0){
dtnu-=tcolmin;
ki=0;
kj=0;
rhu[ki++] = 0.0;
rhv[kj++] = 0.0;
rha[0] = 0.0;
rhb[0] = 0.0;
//******************************//
if(tcol < tcolmin+0.000000001)
{
rc=rbc+icol;
jcol=rc->partnr;
if(jcol==5000000)
{
yycm=rc->ycm-ylc;
zzcm=rc->zcm-zlc;
if(rc->Vycm*yycm >0.0) rc->Vycm=-rc->Vycm;
if(rc->Vzcm*zzcm >0.0) rc->Vzcm=-rc->Vzcm;
}
else
{
if(jcol>70000)
{
printf("ppppppppERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_deltat or tcol");
printf(" rank=%i jcol=%i xsl=%i xer=%i deltaa=%e tcol=%e \n",rank,jcol,xsl,xer,deltaa,tcol);
exit(1);
}
rc2=rbc+jcol;
xxcm=rc->xcm-rc2->xcm;
yycm=rc->ycm-rc2->ycm;
zzcm=rc->zcm-rc2->zcm;
vrx=rc->Vxcm-rc2->Vxcm;
vry=rc->Vycm-rc2->Vycm;
vrz=rc->Vzcm-rc2->Vzcm;
bici=xxcm*vrx+yycm*vry+zzcm*vrz;
bici/=(-4*Rrbc*Rrbc);
rc->Vxcm+= bici*xxcm;
rc->Vycm+= bici*yycm;
rc->Vzcm+= bici*zzcm;
rc2->Vxcm-= bici*xxcm;
rc2->Vycm-= bici*yycm;
rc2->Vzcm-= bici*zzcm;
if(rc->ip == 1)
{
if(rc->xcm <= xsl+Rrbc+0.75)
{
if( (rc2->ip == 0) || ( (rc2->ip == 1) && (rc2->xcm > (xsl+Rrbc+1.0))) )
{
rhu[ki++]=(double)rc->id;
rhu[ki++]=rc->Vxcm;
rhu[ki++]=rc->Vycm;
rhu[ki++]=rc->Vzcm;
}
}
else if(rc->xcm >= xer-Rrbc-0.75)
{
if( (rc2->ip == 0) || ( (rc2->ip == 1) && (rc2->xcm < (xer-Rrbc-1.0))) )
{
rhv[kj++]=(double)rc->id;
rhv[kj++]=rc->Vxcm;
rhv[kj++]=rc->Vycm;
rhv[kj++]=rc->Vzcm;
}
}
}
if(rc2->ip == 1)
{
if(rc2->xcm <= xsl+Rrbc+0.75)
{
if( (rc->ip == 0) || ( (rc->ip == 1) && (rc->xcm > (xsl+Rrbc+1.0))) )
{
rhu[ki++]=(double)rc2->id;
rhu[ki++]=rc2->Vxcm;
rhu[ki++]=rc2->Vycm;
rhu[ki++]=rc2->Vzcm;
}
}
else if(rc2->xcm >= xer-Rrbc-0.75)
{
if( (rc->ip == 0) || ( (rc->ip == 1) && (rc->xcm < (xer-Rrbc-1.0))) )
{
rhv[kj++]=(double)rc2->id;
rhv[kj++]=rc2->Vxcm;
rhv[kj++]=rc2->Vycm;
rhv[kj++]=rc2->Vzcm;
}
}
}
/////////////////////
}
}
////////////////////////////////////////////////////////////////////////
if(ki == 1) ki=0;
if(kj == 1) kj=0;
rhu[0]=(double)ki;
rhv[0]=(double)kj;
rha[0] = 0.0;
rhb[0] = 0.0;
if(numprocs > 1){
MPI_Sendrecv(&rhu[0],ki,MPI_DOUBLE,left,187,&rha[0],5000,MPI_DOUBLE,right,187,comm1d,&status);
MPI_Sendrecv(&rhv[0],kj,MPI_DOUBLE,right,188,&rhb[0],5000,MPI_DOUBLE,left,188,comm1d,&status);
}
ki = (int)rha[0];
kj = 1;
while( kj < ki){
i=(int)rha[kj++];
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip != 1)continue;
if(rc->id == i)
{
rc->Vxcm=rha[kj++];
rc->Vycm=rha[kj++];
rc->Vzcm=rha[kj++];
tcol = tcolmin;
icol=ic;
jcol=6000000;
ic=Nrbcp+10;
}
}
if(ic < Nrbcp+3)
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_ip2");
printf(" rank=%i id=%i xsl=%i xer=%i \n",rank,i,xsl,xer);
exit(1);
}
}
ki = (int)rhb[0];
kj = 1;
while( kj < ki){
i=(int)rhb[kj++];
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip != 1)continue;
if(rc->id == i)
{
rc->Vxcm=rhb[kj++];
rc->Vycm=rhb[kj++];
rc->Vzcm=rhb[kj++];
icol=ic;
jcol=6000000;
tcol = tcolmin;
ic=Nrbcp+10;
}
}
if(ic < Nrbcp+3)
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_ip3");
printf(" rank=%i id=%i xsl=%i xer=%i \n",rank,i,xsl,xer);
exit(1);
}
}
if(tcol < tcolmin+0.000000001)
{
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
if(ic==icol || rc->partnr==icol || ic==jcol || rc->partnr==jcol)
{
rc->coltim=100.0;
rbp=sqrt((ylc-(rc->ycm+rc->Vycm*dtnu))*(ylc-(rc->ycm+rc->Vycm*dtnu))+(zlc-(rc->zcm+rc->Vzcm*dtnu))*(zlc-(rc->zcm+rc->Vzcm*dtnu)));
if(rbp+Rrbc >= Lr0)
{
yycm=rc->ycm-ylc;
zzcm=rc->zcm-zlc;
vry=rc->Vycm;
vrz=rc->Vzcm;
bici=yycm*vry+zzcm*vrz;
if(bici >= 0.0)
{
rsq=yycm*yycm+zzcm*zzcm;
vsq=vry*vry+vrz*vrz;
dic=Lr0-Rrbc;
deltaa=bici*bici-vsq*(rsq-dic*dic);
if(deltaa>0.0)
{
tcol2=(-bici+sqrt(deltaa))/vsq;
if(tcol2 < rc->coltim && tcol2>0.0)
{
rc->coltim=tcol2;
rc->partnr=5000000;
}
}
if(deltaa<0.0 || tcol2<0.0)
{
printf("ffffERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_deltat or tcol2");
printf(" rank=%i ic=%i xsl=%i xer=%i deltaa=%e tcol2=%e \n",rank,ic,xsl,xer,deltaa,tcol2);
exit(1);
}
}
else
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_bici2");
printf(" rank=%i ic=%i xsl=%i xer=%i bici=%e \n",rank,ic,xsl,xer,bici);
exit(1);
}
}
}
}
}
///////////////////////////////////////////////////////////////////
if(tcol < tcolmin+0.000000001)
{
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
if(ic==icol || rc->partnr==icol || ic==jcol || rc->partnr==jcol)
{
for(i=0;i<Nrbcp;i++)
{
rc2=rbc+i;
if(rc2->ip == 2)continue;
if(ic!=i)
{
xxcm=rc->xcm-rc2->xcm;
yycm=rc->ycm-rc2->ycm;
zzcm=rc->zcm-rc2->zcm;
vrx=rc->Vxcm-rc2->Vxcm;
vry=rc->Vycm-rc2->Vycm;
vrz=rc->Vzcm-rc2->Vzcm;
bici=xxcm*vrx+yycm*vry+zzcm*vrz;
if(bici < 0.0)
{
rsq=xxcm*xxcm+yycm*yycm+zzcm*zzcm;
vsq=vrx*vrx+vry*vry+vrz*vrz;
deltaa=bici*bici-vsq*(rsq-4*Rrbc*Rrbc);
if(deltaa>0.0)
{
tcol2=(-bici-sqrt(deltaa))/vsq;
if(tcol2 < rc->coltim && tcol2>0.0)
{
rc->coltim=tcol2;
rc->partnr=i;
}
if(tcol2 < rc2->coltim && tcol2>0.0)
{
rc2->coltim=tcol2;
rc2->partnr=ic;
}
}
}
}
}
}
}
}
/////////////////////////////////////////////////////////
tcol=100.0;
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
if(rc->coltim < tcol)
{
tcol=rc->coltim;
icol=ic;
}
}
MPI_Allreduce(&tcol,&tcolmin,1,MPI_DOUBLE,MPI_MIN,comm1d);
//printf("rankkkkkkkkkkkkkkkkkkkkkkkk22222222222222bbbbbbbbb8=%d tsim=%d tcol=%e dtnu=%e tcolmin=%e \n",rank, tsim, tcol, dtnu, tcolmin);
//////////////////////moving all obstacles with tcolmin ////////////////
if(tcolmin <dtnu)
{
dtnew=tcolmin;
iout=0;
}
else
{
dtnew=dtnu;
iout=1;
}
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
rc->xcm+=rc->Vxcm*dtnew;
rc->ycm+=rc->Vycm*dtnew;
rc->zcm+=rc->Vzcm*dtnew;
rc->coltim-=dtnew;
// if(tsim>495100 && rc->id == 369)
// printf(" rank=%i tsim=%d id=%i ip=%i vxcm=%e vycm=%e vzcm=%e xcm=%e ycm=%e zcm=%e wx=%e wy=%e wz=%e \n",rank,tsim,rc->id,rc->ip,rc->Vxcm,rc->Vycm,rc->Vzcm,rc->xcm,rc->ycm,rc->zcm,rc->wx,rc->wy,rc->wz);
}
///*************************************************************/////
}
//printf("rankkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk2222222222222233333333=%d tcol = %e dt=%e tcolmin=%e iout%d \n",rank, tcol, dt, tcolmin, iout);
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2) continue;
xlo=rc->xcm;
ylo=rc->ycm;
zlo=rc->zcm;
// dtp=dt*0.5;
for(i=0;i<Naf;i++)
{
a=at+i;
if(a->typ == 2)continue;
xl=a->xreal;
yl=a->yreal;
zl=a->zreal;
rbpo=sqrt((xl-xlo)*(xl-xlo)+(yl-ylo)*(yl-ylo)+(zl-zlo)*(zl-zlo));
dtp=dt*0.1;
if(rbpo < Rrbc)
{
rc->Vxcm=-rc->Vxcm;
rc->Vycm=-rc->Vycm;
rc->Vzcm=-rc->Vzcm;
out=0;
j=0;
do
{
rc->xcm+=rc->Vxcm*dtp;
rc->ycm+=rc->Vycm*dtp;
rc->zcm+=rc->Vzcm*dtp;
j++;
xlo=rc->xcm;
ylo=rc->ycm;
zlo=rc->zcm;
rbpo=sqrt((xl-xlo)*(xl-xlo)+(yl-ylo)*(yl-ylo)+(zl-zlo)*(zl-zlo));
rbp=sqrt((ylc-rc->ycm)*(ylc-rc->ycm)+(zlc-rc->zcm)*(zlc-rc->zcm));
if(rbpo >= Rrbc || j>10 || rbp+Rrbc >= Lr0)
{
out=1;
if(rbp+Rrbc >= Lr0)
{
rc->xcm-=rc->Vxcm*dtp;
rc->ycm-=rc->Vycm*dtp;
rc->zcm-=rc->Vzcm*dtp;
}
rc->Vxcm=0.0;
rc->Vycm=0.0;
rc->Vzcm=0.0;
}
}
while(out==0);
//printf("btttttttRRRRRRRRRRROOOOOhjkllgdfOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRa=%d tsim=%d xcm=%e Vxcm=%e wx=%e ip=%d ic=%d id=%d rbpo=%le xer=%d \n",rank,tsim, rc->xcm,rc->Vxcm,rc->wx,rc->ip, ic,rc->id, rbpo,xer);
rbp=sqrt((ylc-rc->ycm)*(ylc-rc->ycm)+(zlc-rc->zcm)*(zlc-rc->zcm));
if( (Rrbc-rbpo) > 0.004 )
{
rc->ip = 2;
// printf("ERRRRRRRRRRRRRRRRROOOOOhjkllgdfOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRa=%d tsim=%d xcm=%e Vxcm=%e wx=%e ip=%d ic=%d id=%d rbpo=%le xer=%d \n",rank,tsim, rc->xcm,rc->Vxcm,rc->wx,rc->ip, ic,rc->id, rbpo,xer);
}
}
}
}
ki=0;
kj=0;
rhu[ki++] = 0.0;
rhv[kj++] = 0.0;
rha[0] = 0.0;
rhb[0] = 0.0;
pi=0;
rpshu[pi++] = 0.0;
pj=0;
rpshv[pj++] = 0.0;
rpsha[0] = 0.0;
rpshb[0] = 0.0;
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2) continue;
if(rc->ip == 0)
{
if(rc->xcm <= xsl+Rrbc+0.5)
{
//printf("rankkkkkkkkkkkkkkkkkkkkkkkkdddddddddddddddddddu=%d tsim=%d xcm=%e dtnu=%e ic=%d id=%d xsl=%d xer=%d \n",rank,tsim, rc->xcm, dtnu, ic,rc->id, xsl,xer);
rc->ip = 1;
rhu[ki++]=(double)rc->id;
rhu[ki++]=rc->xcm;
rhu[ki++]=rc->ycm;
rhu[ki++]=rc->zcm;
rhu[ki++]=rc->Vxcm;
rhu[ki++]=rc->Vycm;
rhu[ki++]=rc->Vzcm;
rhu[ki++]=rc->wx;
rhu[ki++]=rc->wy;
rhu[ki++]=rc->wz;
rhu[ki++]=rc->ex;
rhu[ki++]=rc->ey;
rhu[ki++]=rc->ez;
}
else if(rc->xcm >= xer-Rrbc-0.5)
{
//printf("rankkkkkkkkkkkkkkkkkkkkkkkkdddddddddddddddddddv=%d tsim=%d xcm=%e dtnu=%e ic=%d id=%d xsl=%d xer=%d \n",rank,tsim, rc->xcm, dtnu, ic,rc->id, xsl,xer);
rc->ip = 1;
rhv[kj++]=(double)rc->id;
rhv[kj++]=rc->xcm;
rhv[kj++]=rc->ycm;
rhv[kj++]=rc->zcm;
rhv[kj++]=rc->Vxcm;
rhv[kj++]=rc->Vycm;
rhv[kj++]=rc->Vzcm;
rhv[kj++]=rc->wx;
rhv[kj++]=rc->wy;
rhv[kj++]=rc->wz;
rhv[kj++]=rc->ex;
rhv[kj++]=rc->ey;
rhv[kj++]=rc->ez;
}
}
else if(rc->ip == 1)
{
if((rc->xcm > xsl+Rrbc+0.5) && (rc->xcm < xer-Rrbc-0.5) && (rc->xcm < (xer+xsl)/2.0) )
{
rc->ip = 0;
rpshu[pi++]=(double)rc->id;
// printf("rankkkkkkkkkkkkkkkkkkkkk6666666666666666666666=%d tsim=%d xcm = %e xer=%d xsl=%d ic=%d ip=%d pi=%d id=%d \n",rank,tsim, rc->xcm,xer,xsl,ic,rc->ip,pi,rc->id);
}
if((rc->xcm > xsl+Rrbc+0.5) && (rc->xcm < xer-Rrbc-0.5) && (rc->xcm > (xer+xsl)/2.0))
{
rc->ip = 0;
rpshv[pj++]=(double)rc->id;
// printf("trankkkkkkkkkkkkkkkkkkkkk6666666666666666666666=%d tsim=%d xcm = %e xer=%d xsl=%d ic=%d ip=%d pj=%d id=%d \n",rank,tsim, rc->xcm,xer,xsl,ic,rc->ip,pj,rc->id);
}
}
}
if(ki == 1) ki=0;
if(kj == 1) kj=0;
rhu[0]=(double)ki;
rhv[0]=(double)kj;
rha[0] = 0.0;
rhb[0] = 0.0;
if(pi == 1) pi=0;
if(pj == 1) pj=0;
rpshu[0]=(double)pi;
rpshv[0]=(double)pj;
//printf("rankkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk6666666666666666666666cccccccccccc=%d xer=%d xsl=%d ki=%d kj=%d pi=%d \n",rank, xer,xsl,ki,kj,pi);
if(numprocs > 1){
MPI_Sendrecv(&rhu[0],ki,MPI_DOUBLE,left,189,&rha[0],5000,MPI_DOUBLE,right,189,comm1d,&status);
MPI_Sendrecv(&rhv[0],kj,MPI_DOUBLE,right,190,&rhb[0],5000,MPI_DOUBLE,left,190,comm1d,&status);
MPI_Sendrecv(&rpshu[0],pi,MPI_DOUBLE,left,191,&rpsha[0],5000,MPI_DOUBLE,right,191,comm1d,&status);
MPI_Sendrecv(&rpshv[0],pj,MPI_DOUBLE,right,192,&rpshb[0],5000,MPI_DOUBLE,left,192,comm1d,&status);
}
///*****************************************//
ki = (int)rha[0];
kj = 1;
while( kj < ki){
i=(int)rha[kj++];
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip != 2)continue;
rc->id = i;
rc->ip = 1;
rc->xcm=rha[kj++];
rc->ycm=rha[kj++];
rc->zcm=rha[kj++];
if(rank==numprocs-1) rc->xcm+=(double)L;
rc->Vxcm=rha[kj++];
rc->Vycm=rha[kj++];
rc->Vzcm=rha[kj++];
rc->wx=rha[kj++];
rc->wy=rha[kj++];
rc->wz=rha[kj++];
rc->ex=rha[kj++];
rc->ey=rha[kj++];
rc->ez=rha[kj++];
rc->Momx=0.0;
rc->Momy=0.0;
rc->Momz=0.0;
rc->Amomx=0.0;
rc->Amomy=0.0;
rc->Amomz=0.0;
ic = Nrbcp+10;
}
if(ic < Nrbcp+3)
{
ic = Nrbcp;
Nrbcp++;
rc=rbc+ic;
rc->id = i;
rc->ip = 1;
rc->xcm=rha[kj++];
rc->ycm=rha[kj++];
rc->zcm=rha[kj++];
if(rank==numprocs-1) rc->xcm+=(double)L;
rc->Vxcm=rha[kj++];
rc->Vycm=rha[kj++];
rc->Vzcm=rha[kj++];
rc->wx=rha[kj++];
rc->wy=rha[kj++];
rc->wz=rha[kj++];
rc->ex=rha[kj++];
rc->ey=rha[kj++];
rc->ez=rha[kj++];
rc->Momx=0.0;
rc->Momy=0.0;
rc->Momz=0.0;
rc->Amomx=0.0;
rc->Amomy=0.0;
rc->Amomz=0.0;
}
}
///*****************************************//
ki = (int)rhb[0];
kj = 1;
while( kj < ki){
i=(int)rhb[kj++];
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip != 2)continue;
rc->id = i;
rc->ip = 1;
rc->xcm=rhb[kj++];
rc->ycm=rhb[kj++];
rc->zcm=rhb[kj++];
if(rank==0) rc->xcm-=(double)L;
rc->Vxcm=rhb[kj++];
rc->Vycm=rhb[kj++];
rc->Vzcm=rhb[kj++];
rc->wx=rhb[kj++];
rc->wy=rhb[kj++];
rc->wz=rhb[kj++];
rc->ex=rhb[kj++];
rc->ey=rhb[kj++];
rc->ez=rhb[kj++];
rc->Momx=0.0;
rc->Momy=0.0;
rc->Momz=0.0;
rc->Amomx=0.0;
rc->Amomy=0.0;
rc->Amomz=0.0;
ic = Nrbcp+10;
}
if(ic < Nrbcp+3)
{
ic = Nrbcp;
Nrbcp++;
rc=rbc+ic;
rc->id = i;
rc->ip = 1;
rc->xcm=rhb[kj++];
rc->ycm=rhb[kj++];
rc->zcm=rhb[kj++];
if(rank==0) rc->xcm-=(double)L;
rc->Vxcm=rhb[kj++];
rc->Vycm=rhb[kj++];
rc->Vzcm=rhb[kj++];
rc->wx=rhb[kj++];
rc->wy=rhb[kj++];
rc->wz=rhb[kj++];
rc->ex=rhb[kj++];
rc->ey=rhb[kj++];
rc->ez=rhb[kj++];
rc->Momx=0.0;
rc->Momy=0.0;
rc->Momz=0.0;
rc->Amomx=0.0;
rc->Amomy=0.0;
rc->Amomz=0.0;
}
}
///*****************************************//
ki = (int)rpsha[0];
kj = 1;
while( kj < ki){
i=(int)rpsha[kj++];
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2) continue;
if(rc->id != i)continue;
if(rc->ip != 1)
{
printf("rankkkkkkkkkkkkb=%d xcm=%e xer=%d xsl=%d ic=%d ip=%d kj=%d id=%d \n",rank, rc->xcm,xer,xsl,ic,rc->ip,kj,rc->id);
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_ip2");
printf(" rank=%i ic=%i id=%i ip=%i \n",rank,ic,i,rc->ip);
exit(1);
}
rc->ip = 2;
ic= Nrbcp+10;
}
if(ic < Nrbcp+3)
{
printf("lkeERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_ip5");
printf(" rank=%i id=%i xsl=%i xer=%i \n",rank,i,xsl,xer);
exit(1);
}
}
ki = (int)rpshb[0];
kj = 1;
while( kj < ki){
i=(int)rpshb[kj++];
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2) continue;
if(rc->id != i)continue;
if(rc->ip != 1)
{
printf("rankkkkkkkkkkkkb2=%d xcm=%e xer=%d xsl=%d ic=%d ip=%d kj=%d id=%d \n",rank, rc->xcm,xer,xsl,ic,rc->ip,kj,rc->id);
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_ip43");
printf(" rank=%i ic=%i id=%i ip=%i \n",rank,ic,i,rc->ip);
exit(1);
}
rc->ip = 2;
ic= Nrbcp+10;
}
if(ic < Nrbcp+3)
{
printf("htrERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_ip6");
printf(" rank=%i id=%i xsl=%i xer=%i \n",rank,i,xsl,xer);
exit(1);
}
}
///*****************************************//
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2) continue;
rbp=sqrt((ylc-rc->ycm)*(ylc-rc->ycm)+(zlc-rc->zcm)*(zlc-rc->zcm));
if(rbp+Rrbc >= Lr0)
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_rbp");
printf(" rank=%i ic=%i ip=%i rbp=%e \n",rank,ic,rc->ip,rbp);
exit(1);
}
}
ki=0;
kj=0;
rhu[ki++] = 0.0;
rhv[kj++] = 0.0;
rha[0] = 0.0;
rhb[0] = 0.0;
for(i=0;i<Nrbcp;i++)
{
rc=rbc+i;
if(rc->ip == 2) continue;
if( rc->xcm < xsl-Rrbc-.5)
{
rc->ip = 2;
rhu[ki++]=(double)rc->id;
rhu[ki++]=rc->xcm;
rhu[ki++]=rc->ycm;
rhu[ki++]=rc->zcm;
rhu[ki++]=rc->Vxcm;
rhu[ki++]=rc->Vycm;
rhu[ki++]=rc->Vzcm;
rhu[ki++]=rc->wx;
rhu[ki++]=rc->wy;
rhu[ki++]=rc->wz;
rhu[ki++]=rc->ex;
rhu[ki++]=rc->ey;
rhu[ki++]=rc->ez;
}
else if(rc->xcm > xer+Rrbc+.5)
{
rc->ip = 2;
rhv[kj++]=(double)rc->id;
rhv[kj++]=rc->xcm;
rhv[kj++]=rc->ycm;
rhv[kj++]=rc->zcm;
rhv[kj++]=rc->Vxcm;
rhv[kj++]=rc->Vycm;
rhv[kj++]=rc->Vzcm;
rhv[kj++]=rc->wx;
rhv[kj++]=rc->wy;
rhv[kj++]=rc->wz;
rhv[kj++]=rc->ex;
rhv[kj++]=rc->ey;
rhv[kj++]=rc->ez;
}
}
if(ki == 1) ki=0;
if(kj == 1) kj=0;
rhu[0]=(double)ki;
rhv[0]=(double)kj;
if(numprocs > 1){
MPI_Sendrecv(&rhu[0],ki,MPI_DOUBLE,left,197,&rha[0],5000,MPI_DOUBLE,right,197,comm1d,&status);
MPI_Sendrecv(&rhv[0],kj,MPI_DOUBLE,right,198,&rhb[0],5000,MPI_DOUBLE,left,198,comm1d,&status);
}
///*****************************************//
ki = (int)rha[0];
kj = 1;
while( kj < ki){
i=(int)rha[kj++];
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2) continue;
if(rc->id != i)continue;
rc->ip = 0;
rc->xcm=rha[kj++];
rc->ycm=rha[kj++];
rc->zcm=rha[kj++];
if(rank==numprocs-1) rc->xcm+=(double)L;
rc->Vxcm=rha[kj++];
rc->Vycm=rha[kj++];
rc->Vzcm=rha[kj++];
rc->wx=rha[kj++];
rc->wy=rha[kj++];
rc->wz=rha[kj++];
rc->ex=rha[kj++];
rc->ey=rha[kj++];
rc->ez=rha[kj++];
ic = Nrbcp+10;
}
if(ic < Nrbcp+3)
{
printf("lkeERRRRRfdgRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_ip5");
printf(" rank=%i id=%i xsl=%i xer=%i \n",rank,i,xsl,xer);
exit(1);
}
}
///*****************************************//
ki = (int)rhb[0];
kj = 1;
while( kj < ki){
i=(int)rhb[kj++];
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2) continue;
if(rc->id != i)continue;
rc->ip = 0;
rc->xcm=rhb[kj++];
rc->ycm=rhb[kj++];
rc->zcm=rhb[kj++];
if(rank==0) rc->xcm-=(double)L;
rc->Vxcm=rhb[kj++];
rc->Vycm=rhb[kj++];
rc->Vzcm=rhb[kj++];
rc->wx=rhb[kj++];
rc->wy=rhb[kj++];
rc->wz=rhb[kj++];
rc->ex=rhb[kj++];
rc->ey=rhb[kj++];
rc->ez=rhb[kj++];
ic = Nrbcp+10;
}
if(ic < Nrbcp+3)
{
printf("lkeERRRRRfdgRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_ip5");
printf(" rank=%i id=%i xsl=%i xer=%i \n",rank,i,xsl,xer);
exit(1);
}
}
//***********************************************//
for(i=0;i<Nrbcp;i++)
{
rc=rbc+i;
if(rc->ip == 2) continue;
for(ic=0;ic<Nrbcp;ic++)
{
if(i == ic) continue;
rc2=rbc+ic;
if(rc2->ip == 2) continue;
if(rc->id == rc2->id) continue;
rctoc=sqrt((rc->xcm-rc2->xcm)*(rc->xcm-rc2->xcm)+(rc->ycm-rc2->ycm)*(rc->ycm-rc2->ycm)+(rc->zcm-rc2->zcm)*(rc->zcm-rc2->zcm));
if(rctoc < 2.0*Rrbc)
{
xxcm=rc->xcm-rc2->xcm;
yycm=rc->ycm-rc2->ycm;
zzcm=rc->zcm-rc2->zcm;
if(rc->Vxcm*xxcm < 0.0)rc->Vxcm=-rc->Vxcm;
if(rc->Vycm*yycm < 0.0)rc->Vycm=-rc->Vycm;
if(rc->Vzcm*zzcm < 0.0)rc->Vzcm=-rc->Vzcm;
if(rc2->Vxcm*xxcm > 0.0)rc2->Vxcm=-rc2->Vxcm;
if(rc2->Vycm*yycm > 0.0)rc2->Vycm=-rc2->Vycm;
if(rc2->Vzcm*zzcm > 0.0)rc2->Vzcm=-rc2->Vzcm;
if(rctoc < 1.7*Rrbc)
{
if(rc->id > rc2->id)rc2->ip = 2;
else rc->ip = 2;
}
// exit(1);
}
}
}
}
//***********************************************//
}
double gaussianrandom()
{
double x,y,r2;
do
{
x=-1.+2.*drand48();
y=-1.+2.*drand48();
r2=x*x+y*y;
}
while(r2>=1.0);
return(SRTEMP*y*sqrt(-2.0*log(r2)/r2));
}
void calculateforces()
{
unsigned long int i;
struct aftr *a;
calculateroundspring();
bodyspring();
cenbenspring();
bendbodyforce2();
forcebend11();
for(i=0;i<Naf;i++)
{
a=at+i;
if(a->typ == 0){
a->ax=a->fx;
a->ay=a->fy;
a->az=a->fz;
a->vx+=0.5*dtaf*a->ax;
a->vy+=0.5*dtaf*a->ay;
a->vz+=0.5*dtaf*a->az;
}
a->fx=a->fy=a->fz=0;
}
}
void calculateroundspring()
{
unsigned long int i,j,k;
double dam1x,dam1y,dam1z,d2,d,inv_d,fx,fy,fz,l0;
double pe,tfx,tfy,tfz;
struct aftr *a1,*a2;
tfx=tfy=tfz=0;
for(i=0;i<length;i++)
{
k=i;
j=0;
do
{
if(j==0)
{
a1=(at+k);
a2=(at+k+length*(j+1));
}
else
{
a1=a2;
if(j<npoints-1)
a2=(at+k+length*(j+1));
else
a2=(at+k);
}
j++;
if((a1->typ != 0) && (a2->typ != 0))continue;
if((a2->typ == 2) || (a1->typ == 2))
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRR6");
printf(" rank=%i i=%i a1->typ=%i a2->typ=%i a1->xreal=%le xsl=%i xer=%i \n",rank,i,a1->typ,a2->typ,a1->xreal,xsl,xer);
exit(1);
}
l0=a1->round;
dam1x=(a1->xreal-a2->xreal);
dam1y=(a1->yreal-a2->yreal);
dam1z=(a1->zreal-a2->zreal);
d2=dam1x*dam1x+dam1y*dam1y+dam1z*dam1z;
d=sqrt(d2);
inv_d=1./d;
fx=fy=fz=0;
if(d>0)
{
fx=k0*(d-l0)*dam1x*inv_d*kround;
fy=k0*(d-l0)*dam1y*inv_d*kround;
fz=k0*(d-l0)*dam1z*inv_d*kround;
}
a1->fx+=-fx;
a1->fy+=-fy;
a1->fz+=-fz;
a2->fx+=fx;
a2->fy+=fy;
a2->fz+=fz;
}
while(j<npoints);
}
}
void bodyspring()
{
unsigned long int i,j;
double dam1x,dam1y,dam1z,d2,d,inv_d,fx,fy,fz,l0;
struct aftr *a1,*a2;
i=0;
do
{
for(j=0;j<length-1;j++)
{
a1=at+(i)*length+j;
a2=at+(i)*length+j+1;
if((a1->typ != 0) && (a2->typ != 0))continue;
if((a2->typ == 2) || (a1->typ == 2))
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRR7");
printf(" rank=%i i=%i a1->typ=%i a2->typ=%i a1->xreal=%le xsl=%i xer=%i \n",rank,i,a1->typ,a2->typ,a1->xreal,xsl,xer);
exit(1);
}
l0=a1->body;
dam1x=(a1->xreal-a2->xreal);
dam1y=(a1->yreal-a2->yreal);
dam1z=(a1->zreal-a2->zreal);
//d2=dist2(a1->xreal,a2->xreal,a1->yreal,a2->yreal,a1->zreal,a2->zreal);
d2=dam1x*dam1x+dam1y*dam1y+dam1z*dam1z;
d=sqrt(d2);
inv_d=1./d;
fx=fy=fz=0;
if(d>0)
{
fx=k0*(d-l0)*dam1x*inv_d*kbody;
fy=k0*(d-l0)*dam1y*inv_d*kbody;
fz=k0*(d-l0)*dam1z*inv_d*kbody;
}
a1->fx+=-fx;
a1->fy+=-fy;
a1->fz+=-fz;
a2->fx+=fx;
a2->fy+=fy;
a2->fz+=fz;
}
i++;
}
while(i<npoints);
//printf("\n\n");
//exit(1);
//printf(" %ld\n",num);
}
void cenbenspring()
{
unsigned long int i,j,k,m;
double dam1x,dam1y,dam1z,fx,fy,fz,l0,d2,d,inv_d;
struct aftr *a1,*a2;
for(i=0;i<length;i++)
{
k=i;
//k=length-1;
j=0;
do
{
a1=(at+(k+j*length));
a2=(at+k+length*(j+5));
j++;
if((a1->typ != 0) && (a2->typ != 0))continue;
if((a2->typ == 2) || (a1->typ == 2))
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRR8");
printf(" rank=%i i=%i a1->typ=%i a2->typ=%i a1->xreal=%le xsl=%i xer=%i \n",rank,i,a1->typ,a2->typ,a1->xreal,xsl,xer);
exit(1);
}
l0=a1->cenben;
dam1x=(a1->xreal-a2->xreal);
dam1y=(a1->yreal-a2->yreal);
dam1z=(a1->zreal-a2->zreal);
//d2=dist2(a1->xreal,a2->xreal,a1->yreal,a2->yreal,a1->zreal,a2->zreal);
d2=dam1x*dam1x+dam1y*dam1y+dam1z*dam1z;
d=sqrt(d2);
inv_d=1./d;
fx=fy=fz=0;
if(d>0)
{
fx=k0*(d-l0)*dam1x*inv_d*kcenben;
fy=k0*(d-l0)*dam1y*inv_d*kcenben;
fz=k0*(d-l0)*dam1z*inv_d*kcenben;
}
a1->fx+=-fx;
a1->fy+=-fy;
a1->fz+=-fz;
a2->fx+=fx;
a2->fy+=fy;
a2->fz+=fz;
}
while(j<npoints-5);
//printf("\n\n");
}
//exit(1);
}
void bendbodyforce2()
{
unsigned long int i,j,k;
double mrix,mriy,mriz,mrip1x,mrip1y,mrip1z;
double rix,riy,riz,rip1x,rip1y,rip1z,inv_ri,inv_rip1,inv_mri,inv_mrip1;
double inv_delta,cp,cm,delta,low;
double fx,fy,fz,tfx,tfy,tfz,ttx,tty,ttz,Tx,Ty,Tz;
struct aftr *a1,*a2,*a3;
struct CM *cnmo;
delta=1.e-6;
inv_delta=1./delta;
//cnmo=cm+0;
//printf("%le\n",cnmo->vx);
//exit(1);
j=0;
do
{
low=1.;
for(i=0;i<length-2;i++)
{
a1=at+j*length+i;
a2=at+j*length+i+1;
a3=at+j*length+i+2;
low=(a2->r*a2->r)/(0.78*0.78);
// if(i>9)
// low*=reducedbend;
if((a1->typ != 0) && (a2->typ != 0) && (a3->typ != 0))continue;
if((a1->typ == 2) || (a2->typ == 2) || (a3->typ == 2))
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRR9");
printf(" rank=%i i=%i a1->typ=%i a2->typ=%i a3->typ=%i xsl=%i xer=%i \n",rank,i,a1->typ,a2->typ,a3->typ,xsl,xer);
exit(1);
}
if(a1->typ == 0){
rix=a2->xreal-(a1->xreal+delta);
riy=a2->yreal-a1->yreal;
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal-a2->zreal;
mrix=a2->xreal-(a1->xreal-delta);
mriy=a2->yreal-a1->yreal;
mriz=a2->zreal-a1->zreal;
inv_ri=1./sqrt(rix*rix+riy*riy+riz*riz);
inv_rip1=1./sqrt(rip1x*rip1x+rip1y*rip1y+rip1z*rip1z);
inv_mri=1./sqrt(mrix*mrix+mriy*mriy+mriz*mriz);
rix*=inv_ri;
riy*=inv_ri;
riz*=inv_ri;
rip1x*=inv_rip1;
rip1y*=inv_rip1;
rip1z*=inv_rip1;
mrix*=inv_mri;
mriy*=inv_mri;
mriz*=inv_mri;
mrip1x=rip1x;
mrip1y=rip1y;
mrip1z=rip1z;
cp=rix*rip1x+riy*rip1y+riz*rip1z;
cm=mrix*mrip1x+mriy*mrip1y+mriz*mrip1z;
a1->epx+=0.5*kbend*kmi*(a2->cangle1-cp)*(a2->cangle1-cp)*low;
a1->emx+=0.5*kbend*kmi*(a2->cangle1-cm)*(a2->cangle1-cm)*low;
rix=a2->xreal-a1->xreal;
riy=a2->yreal-a1->yreal-delta;
riz=a2->zreal-a1->zreal;
mrix=a2->xreal-a1->xreal;
mriy=a2->yreal-a1->yreal+delta;
mriz=a2->zreal-a1->zreal;
inv_ri=1./sqrt(rix*rix+riy*riy+riz*riz);
inv_mri=1./sqrt(mrix*mrix+mriy*mriy+mriz*mriz);
rix*=inv_ri;
riy*=inv_ri;
riz*=inv_ri;
mrix*=inv_mri;
mriy*=inv_mri;
mriz*=inv_mri;
cp=rix*rip1x+riy*rip1y+riz*rip1z;
cm=mrix*mrip1x+mriy*mrip1y+mriz*mrip1z;
a1->epy+=0.5*kbend*kmi*(cp-a2->cangle1)*(cp-a2->cangle1)*low;
a1->emy+=0.5*kbend*kmi*(cm-a2->cangle1)*(cm-a2->cangle1)*low;
rix=a2->xreal-a1->xreal;
riy=a2->yreal-a1->yreal;
riz=a2->zreal-a1->zreal-delta;
mrix=a2->xreal-a1->xreal;
mriy=a2->yreal-a1->yreal;
mriz=a2->zreal-a1->zreal+delta;
inv_ri=1./sqrt(rix*rix+riy*riy+riz*riz);
inv_mri=1./sqrt(mrix*mrix+mriy*mriy+mriz*mriz);
rix*=inv_ri;
riy*=inv_ri;
riz*=inv_ri;
mrix*=inv_mri;
mriy*=inv_mri;
mriz*=inv_mri;
cp=rix*rip1x+riy*rip1y+riz*rip1z;
cm=mrix*mrip1x+mriy*mrip1y+mriz*mrip1z;
a1->epz+=0.5*kbend*kmi*(cp-a2->cangle1)*(cp-a2->cangle1)*low;
a1->emz+=0.5*kbend*kmi*(cm-a2->cangle1)*(cm-a2->cangle1)*low;
}
if(a2->typ == 0){
rix=a2->xreal+delta-a1->xreal;
riy=a2->yreal-a1->yreal;
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal-a2->xreal-delta;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal-a2->zreal;
mrix=a2->xreal-delta-a1->xreal;
mriy=a2->yreal-a1->yreal;
mriz=a2->zreal-a1->zreal;
mrip1x=a3->xreal-a2->xreal+delta;
mrip1y=a3->yreal-a2->yreal;
mrip1z=a3->zreal-a2->zreal;
inv_ri=1./sqrt(rix*rix+riy*riy+riz*riz);
inv_rip1=1./sqrt(rip1x*rip1x+rip1y*rip1y+rip1z*rip1z);
inv_mri=1./sqrt(mrix*mrix+mriy*mriy+mriz*mriz);
inv_mrip1=1./sqrt(mrip1x*mrip1x+mrip1y*mrip1y+mrip1z*mrip1z);
rix*=inv_ri;
riy*=inv_ri;
riz*=inv_ri;
rip1x*=inv_rip1;
rip1y*=inv_rip1;
rip1z*=inv_rip1;
mrix*=inv_mri;
mriy*=inv_mri;
mriz*=inv_mri;
mrip1x*=inv_mrip1;
mrip1y*=inv_mrip1;
mrip1z*=inv_mrip1;
cp=rix*rip1x+riy*rip1y+riz*rip1z;
cm=mrix*mrip1x+mriy*mrip1y+mriz*mrip1z;
a2->epx+=0.5*kbend*kmi*(a2->cangle1-cp)*(a2->cangle1-cp)*low;
a2->emx+=0.5*kbend*kmi*(a2->cangle1-cm)*(a2->cangle1-cm)*low;
rix=a2->xreal-a1->xreal;
riy=a2->yreal+delta-a1->yreal;
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal-a2->yreal-delta;
rip1z=a3->zreal-a2->zreal;
mrix=a2->xreal-a1->xreal;
mriy=a2->yreal-delta-a1->yreal;
mriz=a2->zreal-a1->zreal;
mrip1x=a3->xreal-a2->xreal;
mrip1y=a3->yreal-a2->yreal+delta;
mrip1z=a3->zreal-a2->zreal;
inv_ri=1./sqrt(rix*rix+riy*riy+riz*riz);
inv_rip1=1./sqrt(rip1x*rip1x+rip1y*rip1y+rip1z*rip1z);
inv_mri=1./sqrt(mrix*mrix+mriy*mriy+mriz*mriz);
inv_mrip1=1./sqrt(mrip1x*mrip1x+mrip1y*mrip1y+mrip1z*mrip1z);
rix*=inv_ri;
riy*=inv_ri;
riz*=inv_ri;
rip1x*=inv_rip1;
rip1y*=inv_rip1;
rip1z*=inv_rip1;
mrix*=inv_mri;
mriy*=inv_mri;
mriz*=inv_mri;
mrip1x*=inv_mrip1;
mrip1y*=inv_mrip1;
mrip1z*=inv_mrip1;
cp=rix*rip1x+riy*rip1y+riz*rip1z;
cm=mrix*mrip1x+mriy*mrip1y+mriz*mrip1z;
a2->epy+=0.5*kbend*kmi*(cp-a2->cangle1)*(cp-a2->cangle1)*low;
a2->emy+=0.5*kbend*kmi*(cm-a2->cangle1)*(cm-a2->cangle1)*low;
rix=a2->xreal-a1->xreal;
riy=a2->yreal-a1->yreal;
riz=a2->zreal+delta-a1->zreal;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal-a2->zreal-delta;
mrix=a2->xreal-a1->xreal;
mriy=a2->yreal-a1->yreal;
mriz=a2->zreal-delta-a1->zreal;
mrip1x=a3->xreal-a2->xreal;
mrip1y=a3->yreal-a2->yreal;
mrip1z=a3->zreal-a2->zreal+delta;
inv_ri=1./sqrt(rix*rix+riy*riy+riz*riz);
inv_rip1=1./sqrt(rip1x*rip1x+rip1y*rip1y+rip1z*rip1z);
inv_mri=1./sqrt(mrix*mrix+mriy*mriy+mriz*mriz);
inv_mrip1=1./sqrt(mrip1x*mrip1x+mrip1y*mrip1y+mrip1z*mrip1z);
rix*=inv_ri;
riy*=inv_ri;
riz*=inv_ri;
rip1x*=inv_rip1;
rip1y*=inv_rip1;
rip1z*=inv_rip1;
mrix*=inv_mri;
mriy*=inv_mri;
mriz*=inv_mri;
mrip1x*=inv_mrip1;
mrip1y*=inv_mrip1;
mrip1z*=inv_mrip1;
cp=rix*rip1x+riy*rip1y+riz*rip1z;
cm=mrix*mrip1x+mriy*mrip1y+mriz*mrip1z;
a2->epz+=0.5*kbend*kmi*(cp-a2->cangle1)*(cp-a2->cangle1)*low;
a2->emz+=0.5*kbend*kmi*(cm-a2->cangle1)*(cm-a2->cangle1)*low;
}
if(a3->typ == 0){
rix=a2->xreal-a1->xreal;
riy=a2->yreal-a1->yreal;
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal+delta-a2->xreal;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal-a2->zreal;
mrip1x=a3->xreal-delta-a2->xreal;
mrip1y=a3->yreal-a2->yreal;
mrip1z=a3->zreal-a2->zreal;
inv_ri=1./sqrt(rix*rix+riy*riy+riz*riz);
inv_rip1=1./sqrt(rip1x*rip1x+rip1y*rip1y+rip1z*rip1z);
inv_mrip1=1./sqrt(mrip1x*mrip1x+mrip1y*mrip1y+mrip1z*mrip1z);
rix*=inv_ri;
riy*=inv_ri;
riz*=inv_ri;
rip1x*=inv_rip1;
rip1y*=inv_rip1;
rip1z*=inv_rip1;
mrix=rix;
mriy=riy;
mriz=riz;
mrip1x*=inv_mrip1;
mrip1y*=inv_mrip1;
mrip1z*=inv_mrip1;
cp=rix*rip1x+riy*rip1y+riz*rip1z;
cm=mrix*mrip1x+mriy*mrip1y+mriz*mrip1z;
a3->epx+=0.5*kbend*kmi*(a2->cangle1-cp)*(a2->cangle1-cp)*low;
a3->emx+=0.5*kbend*kmi*(a2->cangle1-cm)*(a2->cangle1-cm)*low;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal+delta-a2->yreal;
rip1z=a3->zreal-a2->zreal;
mrip1x=a3->xreal-a2->xreal;
mrip1y=a3->yreal-delta-a2->yreal;
mrip1z=a3->zreal-a2->zreal;
inv_rip1=1./sqrt(rip1x*rip1x+rip1y*rip1y+rip1z*rip1z);
inv_mrip1=1./sqrt(mrip1x*mrip1x+mrip1y*mrip1y+mrip1z*mrip1z);
rip1x*=inv_rip1;
rip1y*=inv_rip1;
rip1z*=inv_rip1;
mrip1x*=inv_mrip1;
mrip1y*=inv_mrip1;
mrip1z*=inv_mrip1;
cp=rix*rip1x+riy*rip1y+riz*rip1z;
cm=mrix*mrip1x+mriy*mrip1y+mriz*mrip1z;
a3->epy+=0.5*kbend*kmi*(cp-a2->cangle1)*(cp-a2->cangle1)*low;
a3->emy+=0.5*kbend*kmi*(cm-a2->cangle1)*(cm-a2->cangle1)*low;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal+delta-a2->zreal;
mrip1x=a3->xreal-a2->xreal;
mrip1y=a3->yreal-a2->yreal;
mrip1z=a3->zreal-delta-a2->zreal;
inv_rip1=1./sqrt(rip1x*rip1x+rip1y*rip1y+rip1z*rip1z);
inv_mrip1=1./sqrt(mrip1x*mrip1x+mrip1y*mrip1y+mrip1z*mrip1z);
rip1x*=inv_rip1;
rip1y*=inv_rip1;
rip1z*=inv_rip1;
mrip1x*=inv_mrip1;
mrip1y*=inv_mrip1;
mrip1z*=inv_mrip1;
cp=rix*rip1x+riy*rip1y+riz*rip1z;
cm=mrix*mrip1x+mriy*mrip1y+mriz*mrip1z;
a3->epz+=0.5*kbend*kmi*(cp-a2->cangle1)*(cp-a2->cangle1)*low;
a3->emz+=0.5*kbend*kmi*(cm-a2->cangle1)*(cm-a2->cangle1)*low;
}
}
//printf("\n\n");
j++;
}
while(j<npoints);
tfx=tfy=tfz=ttx=tty=ttz=Tx=Ty=Tz=0;
for(i=0;i<Naf;i++)
{
a1=at+i;
if(a1->typ == 0){
fx=-(a1->epx-a1->emx)*0.5*inv_delta;
fy=-(a1->epy-a1->emy)*0.5*inv_delta;
fz=-(a1->epz-a1->emz)*0.5*inv_delta;
a1->fx+=fx;
a1->fy+=fy;
a1->fz+=fz;
}
a1->epx=a1->emx=a1->epy=a1->emy=a1->epz=a1->emz=0;
}
//printf("%le\t%le\t%le\n",tfx,tfy,tfz);
//exit(1);
}
void forcebend11()
{
unsigned long int i,j,end,n1,n2,st;
double rix,riy,riz,rip1x,rip1y,rip1z;
double mrix,mriy,mriz,mrip1x,mrip1y,mrip1z;
double rtx,rty,rtz,ux,uy,uz,tx,ty,tz,bx,by,bz,inv_d;
double epx,epy,epz,emx,emy,emz,inv_delta;
double D,g,l0,a,b,delta,fx,tfx,fy,fz,tfy,tfz;
double ttx,tty,ttz,Tx,Ty,Tz,wv11,wv22,h;
double px,py,pz;
struct aftr *a1,*a2,*a3,*s,*a4;
delta=1.e-6;
inv_delta=1./delta;
s=at+1;
//epx=epy=epz=emx=emy=emz=tfx=0;
n1=7*length;
n2=3*length;
end=length-2;
st=2;
//wv11=2.*PI/tlength*1.5;
wv11=wv;
h=gamma1;
a1=at+st;
for(i=st;i<end;i++)
{
if(i==st)
{
a2=a1->helix[0];
a3=a2->helix[1];
}
else
{
a1=a2;
a2=a1->helix[1];
a3=a2->helix[1];
}
if((a1->typ != 0) && (a2->typ != 0) && (a3->typ != 0))continue;
if((a1->typ == 2) || (a2->typ == 2) || (a3->typ == 2))
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRR10");
printf(" rank=%i i=%i a1->typ=%i a2->typ=%i a3->typ=%i xsl=%i xer=%i \n",rank,i,a1->typ,a2->typ,a3->typ,xsl,xer);
exit(1);
}
D=a2->D;
g=avgspcur+h*(sin(wv11*D+fr*totdt));
if(a1->typ == 0){
tx=a1->xreal+delta-a2->xreal;
ty=a1->yreal-a2->yreal;
tz=a1->zreal-a2->zreal;
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
rix=a2->xreal-(a1->xreal+delta);
riy=a2->yreal-a1->yreal;
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal-a2->zreal;
mrix=a2->xreal-(a1->xreal-delta);
mriy=a2->yreal-a1->yreal;
mriz=a2->zreal-a1->zreal;
mrip1x=a3->xreal-a2->xreal;
mrip1y=a3->yreal-a2->yreal;
mrip1z=a3->zreal-a2->zreal;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
a1->epx+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
tx=a1->xreal-delta-a2->xreal;
ty=a1->yreal-a2->yreal;
tz=a1->zreal-a2->zreal;
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
a1->emx+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
if(a2->typ == 0){
rix=a2->xreal+delta-a1->xreal;
riy=a2->yreal-a1->yreal;
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal-(a2->xreal+delta);
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal-a2->zreal;
mrix=a2->xreal-delta-a1->xreal;
mriy=a2->yreal-a1->yreal;
mriz=a2->zreal-a1->zreal;
mrip1x=a3->xreal-a2->xreal+delta;
mrip1y=a3->yreal-a2->yreal;
mrip1z=a3->zreal-a2->zreal;
tx=a1->xreal-(a2->xreal+delta);
ty=a1->yreal-a2->yreal;
tz=a1->zreal-a2->zreal;
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
a2->epx+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
tx=a1->xreal-(a2->xreal-delta);
ty=a1->yreal-a2->yreal;
tz=a1->zreal-a2->zreal;
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
a2->emx+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
if(a3->typ == 0){
rix=a2->xreal-a1->xreal;
riy=a2->yreal-a1->yreal;
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal+delta-a2->xreal;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal-a2->zreal;
mrix=a2->xreal-a1->xreal;
mriy=a2->yreal-a1->yreal;
mriz=a2->zreal-a1->zreal;
mrip1x=a3->xreal-delta-a2->xreal;
mrip1y=a3->yreal-a2->yreal;
mrip1z=a3->zreal-a2->zreal;
tx=a1->xreal-a2->xreal;
ty=a1->yreal-a2->yreal;
tz=a1->zreal-a2->zreal;
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
a3->epx+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
tx=a1->xreal-a2->xreal;
ty=a1->yreal-a2->yreal;
tz=a1->zreal-a2->zreal;
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
a3->emx+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
if(a1->typ == 0){
rix=a2->xreal-a1->xreal;
riy=a2->yreal-(a1->yreal+delta);
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal-a2->zreal;
mrix=a2->xreal-a1->xreal;
mriy=a2->yreal-(a1->yreal-delta);
mriz=a2->zreal-a1->zreal;
mrip1x=a3->xreal-a2->xreal;
mrip1y=a3->yreal-a2->yreal;
mrip1z=a3->zreal-a2->zreal;
tx=a1->xreal-a2->xreal;
ty=a1->yreal+delta-a2->yreal;
tz=a1->zreal-a2->zreal;
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
a1->epy+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
tx=a1->xreal-a2->xreal;
ty=a1->yreal-delta-a2->yreal;
tz=a1->zreal-a2->zreal;
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
a1->emy+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
if(a2->typ == 0){
rix=a2->xreal-a1->xreal;
riy=a2->yreal+delta-a1->yreal;
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal-a2->yreal-delta;
rip1z=a3->zreal-a2->zreal;
mrix=a2->xreal-a1->xreal;
mriy=a2->yreal-delta-a1->yreal;
mriz=a2->zreal-a1->zreal;
mrip1x=a3->xreal-a2->xreal;
mrip1y=a3->yreal-a2->yreal+delta;
mrip1z=a3->zreal-a2->zreal;
tx=a1->xreal-a2->xreal;
ty=a1->yreal-(a2->yreal+delta);
tz=a1->zreal-a2->zreal;
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
a2->epy+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
tx=a1->xreal-a2->xreal;
ty=a1->yreal-(a2->yreal-delta);
tz=a1->zreal-a2->zreal;
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
a2->emy+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
if(a3->typ == 0){
rix=a2->xreal-a1->xreal;
riy=a2->yreal-a1->yreal;
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal+delta-a2->yreal;
rip1z=a3->zreal-a2->zreal;
mrix=a2->xreal-a1->xreal;
mriy=a2->yreal-a1->yreal;
mriz=a2->zreal-a1->zreal;
mrip1x=a3->xreal-a2->xreal;
mrip1y=a3->yreal-delta-a2->yreal;
mrip1z=a3->zreal-a2->zreal;
tx=a1->xreal-a2->xreal;
ty=a1->yreal-a2->yreal;
tz=a1->zreal-a2->zreal;
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
a3->epy+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
a3->emy+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
if(a1->typ == 0){
rix=a2->xreal-a1->xreal;
riy=a2->yreal-a1->yreal;
riz=a2->zreal-(a1->zreal+delta);
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal-a2->zreal;
mrix=a2->xreal-a1->xreal;
mriy=a2->yreal-a1->yreal;
mriz=a2->zreal-(a1->zreal-delta);
mrip1x=a3->xreal-a2->xreal;
mrip1y=a3->yreal-a2->yreal;
mrip1z=a3->zreal-a2->zreal;
tx=a1->xreal-a2->xreal;
ty=a1->yreal-a2->yreal;
tz=a1->zreal+delta-a2->zreal;
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
a1->epz+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
tx=a1->xreal-a2->xreal;
ty=a1->yreal-a2->yreal;
tz=a1->zreal-delta-a2->zreal;
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
a1->emz+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
if(a2->typ == 0){
rix=a2->xreal-a1->xreal;
riy=a2->yreal-a1->yreal;
riz=a2->zreal+delta-a1->zreal;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal-a2->zreal-delta;
mrix=a2->xreal-a1->xreal;
mriy=a2->yreal-a1->yreal;
mriz=a2->zreal-delta-a1->zreal;
mrip1x=a3->xreal-a2->xreal;
mrip1y=a3->yreal-a2->yreal;
mrip1z=a3->zreal-a2->zreal+delta;
tx=a1->xreal-a2->xreal;
ty=a1->yreal-a2->yreal;
tz=a1->zreal-(a2->zreal+delta);
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
a2->epz+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
tx=a1->xreal-a2->xreal;
ty=a1->yreal-a2->yreal;
tz=a1->zreal-(a2->zreal-delta);
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
a2->emz+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
if(a3->typ == 0){
rix=a2->xreal-a1->xreal;
riy=a2->yreal-a1->yreal;
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal+delta-a2->zreal;
mrix=a2->xreal-a1->xreal;
mriy=a2->yreal-a1->yreal;
mriz=a2->zreal-a1->zreal;
mrip1x=a3->xreal-a2->xreal;
mrip1y=a3->yreal-a2->yreal;
mrip1z=a3->zreal-delta-a2->zreal;
tx=a1->xreal-a2->xreal;
ty=a1->yreal-a2->yreal;
tz=a1->zreal-a2->zreal;
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
a3->epz+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
a3->emz+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
//printf("%ld\t%ld\t%ld\t%ld\t%ld\n",a1->n,a2->n,a3->n,(at+a2->n1)->n,(at+a2->n2)->n);
//printf("%le\t%le\t%le\t%le\t%le\t%le\n",a2->xreal,a2->yreal,a2->zreal,ux,uy,uz);
}
//exit(1);
//D=0;
a1=at+st;
for(i=st;i<end;i++)
{
if(i==st)
{
a2=a1->helix[0];
a3=a2->helix[1];
}
else
{
a1=a2;
a2=a1->helix[1];
a3=a2->helix[1];
}
if(((at+a2->n1)->typ != 0) && ((at+a2->n2)->typ != 0) )continue;
if(((at+a2->n1)->typ == 2) || ((at+a2->n2)->typ == 2) || (a1->typ == 2) || (a2->typ == 2) || (a3->typ == 2))
{
printf("ERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRR12");
printf(" rank=%i i=%i a1->typ=%i a2->typ=%i a3->typ=%i (at+a2->n1)->typ=%i (at+a2->n2)->typ=%i \n",rank,i,a1->typ,a2->typ,a3->typ,(at+a2->n1)->typ,(at+a2->n2)->typ);
exit(1);
}
D=a2->D;
h=gamma1;
g=avgspcur+h*(sin(wv11*D+fr*totdt));
rix=a2->xreal-a1->xreal;
riy=a2->yreal-a1->yreal;
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal-a2->zreal;
mrix=a2->xreal-a1->xreal;
mriy=a2->yreal-a1->yreal;
mriz=a2->zreal-a1->zreal;
mrip1x=a3->xreal-a2->xreal;
mrip1y=a3->yreal-a2->yreal;
mrip1z=a3->zreal-a2->zreal;
tx=a1->xreal-a2->xreal;
ty=a1->yreal-a2->yreal;
tz=a1->zreal-a2->zreal;
l0=sqrt(tx*tx+ty*ty+tz*tz);
a=cos(g*l0);
b=sin(g*l0);
if((at+a2->n1)->typ == 0){
bx=(at+a2->n1)->xreal+delta-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
(at+a2->n1)->epx+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
bx=(at+a2->n1)->xreal-delta-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
(at+a2->n1)->emx+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
if((at+a2->n2)->typ == 0){
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal-delta;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
(at+a2->n2)->epx+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal+delta;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
(at+a2->n2)->emx+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
rix=a2->xreal-a1->xreal;
riy=a2->yreal-a1->yreal;
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal-a2->zreal;
mrix=a2->xreal-a1->xreal;
mriy=a2->yreal-a1->yreal;
mriz=a2->zreal-a1->zreal;
mrip1x=a3->xreal-a2->xreal;
mrip1y=a3->yreal-a2->yreal;
mrip1z=a3->zreal-a2->zreal;
if((at+a2->n1)->typ == 0){
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal+delta-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
(at+a2->n1)->epy+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-delta-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
(at+a2->n1)->emy+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
if((at+a2->n2)->typ == 0){
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal-delta;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
(at+a2->n2)->epy+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal+delta;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
(at+a2->n2)->emy+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
rix=a2->xreal-a1->xreal;
riy=a2->yreal-a1->yreal;
riz=a2->zreal-a1->zreal;
rip1x=a3->xreal-a2->xreal;
rip1y=a3->yreal-a2->yreal;
rip1z=a3->zreal-a2->zreal;
mrix=a2->xreal-a1->xreal;
mriy=a2->yreal-a1->yreal;
mriz=a2->zreal-a1->zreal;
mrip1x=a3->xreal-a2->xreal;
mrip1y=a3->yreal-a2->yreal;
mrip1z=a3->zreal-a2->zreal;
if((at+a2->n1)->typ == 0){
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal+delta-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
(at+a2->n1)->epz+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-delta-(at+a2->n2)->zreal;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
(at+a2->n1)->emz+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
if((at+a2->n2)->typ == 0){
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal-delta;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
rtx=rip1x-(u11(ux,uy,uz,a,b)*rix+u12(ux,uy,uz,a,b)*riy+u13(ux,uy,uz,a,b)*riz);
rty=rip1y-(u21(ux,uy,uz,a,b)*rix+u22(ux,uy,uz,a,b)*riy+u23(ux,uy,uz,a,b)*riz);
rtz=rip1z-(u31(ux,uy,uz,a,b)*rix+u32(ux,uy,uz,a,b)*riy+u33(ux,uy,uz,a,b)*riz);
(at+a2->n2)->epz+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
bx=(at+a2->n1)->xreal-(at+a2->n2)->xreal;
by=(at+a2->n1)->yreal-(at+a2->n2)->yreal;
bz=(at+a2->n1)->zreal-(at+a2->n2)->zreal+delta;
ux=crossx(tx,ty,tz,bx,by,bz);
uy=crossy(tx,ty,tz,bx,by,bz);
uz=crossz(tx,ty,tz,bx,by,bz);
inv_d=1./sqrt(ux*ux+uy*uy+uz*uz);
ux*=inv_d;
uy*=inv_d;
uz*=inv_d;
rtx=mrip1x-(u11(ux,uy,uz,a,b)*mrix+u12(ux,uy,uz,a,b)*mriy+u13(ux,uy,uz,a,b)*mriz);
rty=mrip1y-(u21(ux,uy,uz,a,b)*mrix+u22(ux,uy,uz,a,b)*mriy+u23(ux,uy,uz,a,b)*mriz);
rtz=mrip1z-(u31(ux,uy,uz,a,b)*mrix+u32(ux,uy,uz,a,b)*mriy+u33(ux,uy,uz,a,b)*mriz);
(at+a2->n2)->emz+=0.5*kbend*knew*(rtx*rtx+rty*rty+rtz*rtz);
}
}
a1=at+st;
tfx=tfy=tfz=Tx=Ty=Tz=0;
for(i=st;i<end;i++)
{
if(i==st)
{
a2=a1->helix[0];
a3=a2->helix[1];
}
else
{
a1=a2;
a2=a1->helix[1];
a3=a2->helix[1];
}
if((at+a2->n1)->typ == 0){
fx=-((at+a2->n1)->epx-(at+a2->n1)->emx)*0.5*inv_delta;
fy=-((at+a2->n1)->epy-(at+a2->n1)->emy)*0.5*inv_delta;
fz=-((at+a2->n1)->epz-(at+a2->n1)->emz)*0.5*inv_delta;
(at+a2->n1)->fx+=fx;
(at+a2->n1)->fy+=fy;
(at+a2->n1)->fz+=fz;
tfx+=fx;
tfy+=fy;
tfz+=fz;
ttx=crossx((at+a2->n1)->xreal,(at+a2->n1)->yreal,(at+a2->n1)->zreal,fx,fy,fz);
tty=crossy((at+a2->n1)->xreal,(at+a2->n1)->yreal,(at+a2->n1)->zreal,fx,fy,fz);
ttz=crossz((at+a2->n1)->xreal,(at+a2->n1)->yreal,(at+a2->n1)->zreal,fx,fy,fz);
Tx+=ttx;
Ty+=tty;
Tz+=ttz;
}
(at+a2->n1)->epx=(at+a2->n1)->emx=(at+a2->n1)->epy=(at+a2->n1)->emy=(at+a2->n1)->epz=(at+a2->n1)->emz=0;
if((at+a2->n2)->typ == 0){
fx=-((at+a2->n2)->epx-(at+a2->n2)->emx)*0.5*inv_delta;
fy=-((at+a2->n2)->epy-(at+a2->n2)->emy)*0.5*inv_delta;
fz=-((at+a2->n2)->epz-(at+a2->n2)->emz)*0.5*inv_delta;
(at+a2->n2)->fx+=fx;
(at+a2->n2)->fy+=fy;
(at+a2->n2)->fz+=fz;
tfx+=fx;
tfy+=fy;
tfz+=fz;
ttx=crossx((at+a2->n2)->xreal,(at+a2->n2)->yreal,(at+a2->n2)->zreal,fx,fy,fz);
tty=crossy((at+a2->n2)->xreal,(at+a2->n2)->yreal,(at+a2->n2)->zreal,fx,fy,fz);
ttz=crossz((at+a2->n2)->xreal,(at+a2->n2)->yreal,(at+a2->n2)->zreal,fx,fy,fz);
Tx+=ttx;
Ty+=tty;
Tz+=ttz;
}
(at+a2->n2)->epx=(at+a2->n2)->emx=(at+a2->n2)->epy=(at+a2->n2)->emy=(at+a2->n2)->epz=(at+a2->n2)->emz=0;
}
a2=at+st;
j=1;
for(i=st;i<end+2;i++)
{
a1=a2;
if(a1->typ == 0){
fx=-(a1->epx-a1->emx)*0.5*inv_delta;
fy=-(a1->epy-a1->emy)*0.5*inv_delta;
fz=-(a1->epz-a1->emz)*0.5*inv_delta;
a1->fx+=fx;
a1->fy+=fy;
a1->fz+=fz;
tfx+=fx;
tfy+=fy;
tfz+=fz;
ttx=crossx(a1->xreal,a1->yreal,a1->zreal,fx,fy,fz);
tty=crossy(a1->xreal,a1->yreal,a1->zreal,fx,fy,fz);
ttz=crossz(a1->xreal,a1->yreal,a1->zreal,fx,fy,fz);
Tx+=ttx;
Ty+=tty;
Tz+=ttz;
}
a1->epx=a1->emx=a1->epy=a1->emy=a1->epz=a1->emz=0;
if(i==st)
a2=a1->helix[0];
if(i>st && i<end+1)
a2=a1->helix[1];
}
}
void collisionAT(void)
{
unsigned long int i,n,j,ki=0,kj=0,Nsolpa,ic;
double xl,yl,zl,inv_j;
double vxs,vys,vzs,v2;
double u11,u12,u13,u21,u22,u23,u31,u32,u33,d,inv_d;
double shift1,shift2,shift3,tr,px,py,pz,rpc,rbp;
struct sol *srd;
struct aftr *a;
struct wallp *pwl;
struct redcell *rc;
static double *hux;
static double *hvx;
static double *hax;
static double *hbx;
static int n2t = 0;
if(n2t == 0){
hux = (double *) malloc(sizeof(double)*(600000));
if(hux == NULL){
printf("\n err_message(2)");
exit(1);
}
hvx = (double *) malloc(sizeof(double)*(600000));
if(hvx == NULL) exit(1);
hax = (double *) malloc(sizeof(double)*(600000));
if(hax == NULL) exit(1);
hbx = (double *) malloc(sizeof(double)*(600000));
if(hbx == NULL) exit(1);
n2t++;
}
hux[ki++] = 0.0;
hvx[0] = 0.0;
hvx[1] = 0.0;
hbx[0] = 0.0;
hbx[1] = 0.0;
hax[0] = 0.0;
Nsolpa = Nsolp;
wallpos();
for(i=0;i<Nsolp;i++)
{
srd=ss+i;
srd->pn=1000000000;
}
for(i=0;i<Naf;i++)
{
a=at+i;
a->pn=1000000000;
}
for(i=0;i<Nsolw;i++)
{
pwl=ww+i;
pwl->pn=1000000000;
}
for(i=0;i<Lr1*Lr1*L/numprocs;i++)
{
N_VCM[i].Vx=N_VCM[i].Vy=N_VCM[i].Vz=0;
N_VCM[i].Vnx=N_VCM[i].Vny=N_VCM[i].Vnz=0;
N_VCM[i].xcm=N_VCM[i].ycm=N_VCM[i].zcm=0;
N_VCM[i].xx=N_VCM[i].yy=N_VCM[i].zz=0;
N_monocell[i]=0;
N_VCM[i].S=1.;
N_VCM[i].r11=N_VCM[i].r12=N_VCM[i].r13=0;
N_VCM[i].r21=N_VCM[i].r22=N_VCM[i].r23=0;
N_VCM[i].r31=N_VCM[i].r32=N_VCM[i].r33=0;
N_VCM[i].crx=N_VCM[i].cry=N_VCM[i].crz=0;
}
if(rank==0){
shift2=(drand48()-0.5)*(double)cellsize;
shift3=(drand48()-0.5)*(double)cellsize;
shift1=(drand48()-0.5)*(double)cellsize;
}
if(numprocs > 1)
{
MPI_Bcast(&shift1,1,MPI_DOUBLE,0,comm1d);
MPI_Bcast(&shift2,1,MPI_DOUBLE,0,comm1d);
MPI_Bcast(&shift3,1,MPI_DOUBLE,0,comm1d);
}
//////////////////////////////
for(i=0;i<Nsolp;i++)
{
srd=ss+i;
xl=srd->x;
if(xl < -1000) continue;
if(xl < xsl) printf("\n rankkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk=%i i=%i xl=%e Nslop=%i \n",rank,i,xl,Nsolp);
if(xl >= xer) printf("\n rankkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk22222222222222=%i i=%i xl=%e Nsolp=%i \n",rank,i,xl,Nsolp);
xl=srd->x-shift1;
yl=srd->y-shift2;
zl=srd->z-shift3;
/* if(xl<0)
xl+=(double)L;
if(xl>=(double)L)
xl-=(double)L; */
/*
if(yl<0)
yl+=(double)Lr1;
if(yl>=(double)Lr1)
yl-=(double)Lr1;
if(zl<0)
zl+=(double)Lr1;
if(zl>=(double)Lr1)
zl-=(double)Lr1;
*/
if(yl<0 || yl>=(double)Lr1)
{
printf("\n err_message(300)");//////////
exit(1);
}
if(zl<0 || zl>=(double)Lr1)
{
printf("\n err_message(301)");/////////
exit(1);
}
srd->xa=xl;
srd->ya=yl;
srd->za=zl;
if(shift1 > 0)
{
if(xl < xsl)
{
hux[ki++]=(double)i;
hux[ki++]=xl;
hux[ki++]=yl;
hux[ki++]=zl;
hux[ki++]= srd->vx;
hux[ki++]= srd->vy;
hux[ki++]= srd->vz;
}
else if(xl >= xer)
{
printf("\n err_message(3)");
exit(1);
}
}
else if(shift1 <= 0)
{
if(xl >= xer)
{
hux[ki++]=(double)i;
hux[ki++]=xl;
hux[ki++]=yl;
hux[ki++]=zl;
hux[ki++]= srd->vx;
hux[ki++]= srd->vy;
hux[ki++]= srd->vz;
}
else if(xl < xsl)
{
printf("\n err_message(4)");///////
exit(1);
}
}
}
if(ki == 1) ki=0;
hux[0]=(double)ki;
if(numprocs > 1){
if(shift1 > 0)
MPI_Sendrecv(&hux[0],ki,MPI_DOUBLE,left,7,&hax[0],600000,MPI_DOUBLE,right,7,comm1d,&status);
else
MPI_Sendrecv(&hux[0],ki,MPI_DOUBLE,right,8,&hax[0],600000,MPI_DOUBLE,left,8,comm1d,&status);
// MPI_Sendrecv(&hvp[0],kj,MPI_DOUBLE,right,4,&hbp[0],600000,MPI_DOUBLE,left,4,comm1d,&status);
}
ki = (int)hax[0];
kj = 1;
while( kj < ki){
j = Nsolpa ;
Nsolpa++;
srd=ss+j;
srd->n=(int)hax[kj++];
srd->xa=hax[kj++];
if(shift1 > 0 && rank==numprocs-1) srd->xa+=(double)L;
if(shift1 < 0 && rank==0) srd->xa-=(double)L;
srd->ya=hax[kj++];
srd->za=hax[kj++];
srd->vx=hax[kj++];
srd->vy=hax[kj++];
srd->vz=hax[kj++];
}
for(i=0;i<Nsolpa;i++)
{
srd=ss+i;
xl=srd->x;
if(i<Nsolp && xl < -1000) continue;
xl=srd->xa;
if(xl < xsl || xl >= xer) continue;
yl=srd->ya;
zl=srd->za;
n=(unsigned long int)(floor(xl-xsl)+(double)(L/numprocs)*(floor(yl)+(double)Lr1*floor(zl)));
if(n > Lr1*Lr1*L/numprocs) printf("\n afrankkknnnnnnnnnnnnnpppppppppppppppppppp=%i n=%i xl=%e xsl=%i xer=%i \n",rank,n,xl,xsl,xer);
srd->pn=n;
N_monocell[n]+=1;
N_VCM[n].Vx+=srd->vx;
N_VCM[n].Vy+=srd->vy;
N_VCM[n].Vz+=srd->vz;
N_VCM[n].xcm+=srd->xa;
N_VCM[n].ycm+=srd->ya;
N_VCM[n].zcm+=srd->za;
}
tbf=0;
for(i=0;i<Naf;i++)
{
a=at+i;
if(a->typ ==2 )continue;
xl=a->xreal-shift1;
//xl=a->x;
if(xl < xsl || xl >= xer)continue;
yl=a->yreal-shift2;
zl=a->zreal-shift3;
a->vx+=vaddx;
a->vy+=vaddy;
a->vz+=vaddz;
tbf+=a->vx*a->vx+a->vy*a->vy+a->vz*a->vz;
/* if(xl<0)
xl+=(double)L;
if(xl>=(double)L)
xl-=(double)L; */
/* if(yl<0)
yl+=(double)Lr1;
if(yl>=(double)Lr1)
yl-=(double)Lr1;
if(zl<0)
zl+=(double)Lr1;
if(zl>=(double)Lr1)
zl-=(double)Lr1;
*/
if(yl<0 || yl>=(double)Lr1)
{
printf("\n err_message(302)");
exit(1);
}
if(zl<0 || zl>=(double)Lr1)
{
printf("\n err_message(303)");
exit(1);
}
n=(unsigned long int)(floor(xl-xsl)+(double)(L/numprocs)*(floor(yl)+(double)Lr1*floor(zl)));
if(n > Lr1*Lr1*L/numprocs) printf("\n afrankkknnnnnnnnnnnnnnnnnnnnnnnnnnnn=%i n=%i xl=%e xsl=%i xer=%i \n",rank,n,xl,xsl,xer);
a->pn=n;
N_monocell[n]+=1;
a->xa=xl;
a->ya=yl;
a->za=zl;
N_VCM[n].Vx+=a->vx;
N_VCM[n].Vy+=a->vy;
N_VCM[n].Vz+=a->vz;
N_VCM[n].xcm+=a->xa;
N_VCM[n].ycm+=a->ya;
N_VCM[n].zcm+=a->za;
}
tbf*=inv_Naf;
for(i=0;i<Nsolw;i++)
{
pwl=ww+i;
xl=pwl->x-shift1;;
yl=pwl->y-shift2;
zl=pwl->z-shift3;
pwl->xa=xl;
pwl->ya=yl;
pwl->za=zl;
if(yl<0 || yl>=(double)Lr1) continue;
if(zl<0 || zl>=(double)Lr1) continue;
if(xl < xsl || xl >= xer) continue;
n=(unsigned long int)(floor(xl-xsl)+(double)(L/numprocs)*(floor(yl)+(double)Lr1*floor(zl)));
if(n > Lr1*Lr1*L/numprocs) printf("\n afrankkknnnnnnnnnnnnnppppppppppppppppppppwwwall=%i n=%i xl=%e xsl=%i xer=%i \n",rank,n,xl,xsl,xer);
pwl->pn=n;
N_monocell[n]+=1;
N_VCM[n].Vx+=pwl->vx;
N_VCM[n].Vy+=pwl->vy;
N_VCM[n].Vz+=pwl->vz;
N_VCM[n].xcm+=pwl->xa;
N_VCM[n].ycm+=pwl->ya;
N_VCM[n].zcm+=pwl->za;
}
for(i=0;i<Lr1*Lr1*L/numprocs;i++)
{
j=N_monocell[i];
if(j>0)
{
inv_j=1./(double)j;
N_VCM[i].xcm*=inv_j;
N_VCM[i].ycm*=inv_j;
N_VCM[i].zcm*=inv_j;
N_VCM[i].Vx*=inv_j;
N_VCM[i].Vy*=inv_j;
N_VCM[i].Vz*=inv_j;
N_VCM[i].S=inv_j;
}
}
for(i=0;i<Nsolpa;i++)
{
srd=ss+i;
n=srd->pn;
if(n > 10000000)continue;
xl=srd->x;
if(i<Nsolp && xl < -1000)
{
printf("\n err_message(5)");
exit(1);
}
xl=srd->xa;
if(xl < xsl || xl >= xer)
{
printf("\n rankkknnnnnnnnnnnnnnnnnnnnnnnnnnnn=%i n=%i i=%i xl=%e xsl=%i xer=%i \n",rank,n,i,xl,xsl,xer);
printf("\n err_messagee(6)");
exit(1);
}
vxs=gaussianrandom();
vys=gaussianrandom();
vzs=gaussianrandom();
N_VCM[n].Vnx+=vxs*N_VCM[n].S;
N_VCM[n].Vny+=vys*N_VCM[n].S;
N_VCM[n].Vnz+=vzs*N_VCM[n].S;
xl=srd->xa-N_VCM[n].xcm;
yl=srd->ya-N_VCM[n].ycm;
zl=srd->za-N_VCM[n].zcm;
N_VCM[n].crx+=crossx(xl,yl,zl,(srd->vx-vxs),(srd->vy-vys),(srd->vz-vzs));
N_VCM[n].cry+=crossy(xl,yl,zl,(srd->vx-vxs),(srd->vy-vys),(srd->vz-vzs));
N_VCM[n].crz+=crossz(xl,yl,zl,(srd->vx-vxs),(srd->vy-vys),(srd->vz-vzs));
N_VCM[n].r11+=yl*yl+zl*zl;
N_VCM[n].r12-=xl*yl;
N_VCM[n].r13-=xl*zl;
N_VCM[n].r21-=xl*yl;
N_VCM[n].r22+=xl*xl+zl*zl;
N_VCM[n].r23-=yl*zl;
N_VCM[n].r31-=xl*zl;
N_VCM[n].r32-=yl*zl;
N_VCM[n].r33+=xl*xl+yl*yl;
srd->vx=vxs;
srd->vy=vys;
srd->vz=vzs;
}
for(i=0;i<Naf;i++)
{
a=at+i;
n=a->pn;
if(n > 10000000) continue;
//if(a->typ != 0 && a->typ != 1)continue;
//if( a->xa < xsl || a->xa >= xer)continue;
vxs=gaussianrandom();
vys=gaussianrandom();
vzs=gaussianrandom();
N_VCM[n].Vnx+=vxs*N_VCM[n].S;
N_VCM[n].Vny+=vys*N_VCM[n].S;
N_VCM[n].Vnz+=vzs*N_VCM[n].S;
xl=a->xa-N_VCM[n].xcm;
yl=a->ya-N_VCM[n].ycm;
zl=a->za-N_VCM[n].zcm;
N_VCM[n].crx+=crossx(xl,yl,zl,(a->vx-vxs),(a->vy-vys),(a->vz-vzs));
N_VCM[n].cry+=crossy(xl,yl,zl,(a->vx-vxs),(a->vy-vys),(a->vz-vzs));
N_VCM[n].crz+=crossz(xl,yl,zl,(a->vx-vxs),(a->vy-vys),(a->vz-vzs));
N_VCM[n].r11+=yl*yl+zl*zl;
N_VCM[n].r12-=xl*yl;
N_VCM[n].r13-=xl*zl;
N_VCM[n].r21-=xl*yl;
N_VCM[n].r22+=xl*xl+zl*zl;
N_VCM[n].r23-=yl*zl;
N_VCM[n].r31-=xl*zl;
N_VCM[n].r32-=yl*zl;
N_VCM[n].r33+=xl*xl+yl*yl;
a->vx=vxs;
a->vy=vys;
a->vz=vzs;
}
for(i=0;i<Nsolw;i++)
{
pwl=ww+i;
n=pwl->pn;
if(n > 10000000)continue;
vxs=gaussianrandom();
vys=gaussianrandom();
vzs=gaussianrandom();
N_VCM[n].Vnx+=vxs*N_VCM[n].S;
N_VCM[n].Vny+=vys*N_VCM[n].S;
N_VCM[n].Vnz+=vzs*N_VCM[n].S;
xl=pwl->xa-N_VCM[n].xcm;
yl=pwl->ya-N_VCM[n].ycm;
zl=pwl->za-N_VCM[n].zcm;
N_VCM[n].crx+=crossx(xl,yl,zl,(pwl->vx-vxs),(pwl->vy-vys),(pwl->vz-vzs));
N_VCM[n].cry+=crossy(xl,yl,zl,(pwl->vx-vxs),(pwl->vy-vys),(pwl->vz-vzs));
N_VCM[n].crz+=crossz(xl,yl,zl,(pwl->vx-vxs),(pwl->vy-vys),(pwl->vz-vzs));
N_VCM[n].r11+=yl*yl+zl*zl;
N_VCM[n].r12-=xl*yl;
N_VCM[n].r13-=xl*zl;
N_VCM[n].r21-=xl*yl;
N_VCM[n].r22+=xl*xl+zl*zl;
N_VCM[n].r23-=yl*zl;
N_VCM[n].r31-=xl*zl;
N_VCM[n].r32-=yl*zl;
N_VCM[n].r33+=xl*xl+yl*yl;
pwl->vrx=vxs;
pwl->vry=vys;
pwl->vrz=vzs;
}
for(n=0;n<Lr1*Lr1*L/numprocs;n++)
{
d=0;
d+=N_VCM[n].r11*(N_VCM[n].r22*N_VCM[n].r33-N_VCM[n].r23*N_VCM[n].r32);
d-=N_VCM[n].r12*(N_VCM[n].r21*N_VCM[n].r33-N_VCM[n].r23*N_VCM[n].r31);
d+=N_VCM[n].r13*(N_VCM[n].r21*N_VCM[n].r32-N_VCM[n].r22*N_VCM[n].r31);
if(d*d>1.e-6 && N_monocell[n]>2)
{
u11=(N_VCM[n].r22*N_VCM[n].r33)-(N_VCM[n].r23*N_VCM[n].r32);
u12=(N_VCM[n].r13*N_VCM[n].r32)-(N_VCM[n].r12*N_VCM[n].r33);
u13=(N_VCM[n].r12*N_VCM[n].r23)-(N_VCM[n].r13*N_VCM[n].r22);
u21=(N_VCM[n].r23*N_VCM[n].r31)-(N_VCM[n].r21*N_VCM[n].r33);
u22=(N_VCM[n].r11*N_VCM[n].r33)-(N_VCM[n].r13*N_VCM[n].r31);
u23=(N_VCM[n].r13*N_VCM[n].r21)-(N_VCM[n].r11*N_VCM[n].r23);
u31=(N_VCM[n].r21*N_VCM[n].r32)-(N_VCM[n].r22*N_VCM[n].r31);
u32=(N_VCM[n].r12*N_VCM[n].r31)-(N_VCM[n].r11*N_VCM[n].r32);
u33=(N_VCM[n].r11*N_VCM[n].r22)-(N_VCM[n].r12*N_VCM[n].r21);
inv_d=1./d;
N_VCM[n].r11=u11*inv_d;
N_VCM[n].r12=u12*inv_d;
N_VCM[n].r13=u13*inv_d;
N_VCM[n].r21=u21*inv_d;
N_VCM[n].r22=u22*inv_d;
N_VCM[n].r23=u23*inv_d;
N_VCM[n].r31=u31*inv_d;
N_VCM[n].r32=u32*inv_d;
N_VCM[n].r33=u33*inv_d;
}
else
{
d=(N_VCM[n].r11+N_VCM[n].r22+N_VCM[n].r33)*(N_VCM[n].r11+N_VCM[n].r22+N_VCM[n].r33);
if(d*d>1.e-6 && N_monocell[n]==2)
{
tr=4./d;
N_VCM[n].r11*=tr;
N_VCM[n].r12*=tr;
N_VCM[n].r13*=tr;
N_VCM[n].r21*=tr;
N_VCM[n].r22*=tr;
N_VCM[n].r23*=tr;
N_VCM[n].r31*=tr;
N_VCM[n].r32*=tr;
N_VCM[n].r33*=tr;
}
else
{
N_VCM[n].r11=0;
N_VCM[n].r12=0;
N_VCM[n].r13=0;
N_VCM[n].r21=0;
N_VCM[n].r22=0;
N_VCM[n].r23=0;
N_VCM[n].r31=0;
N_VCM[n].r32=0;
N_VCM[n].r33=0;
}
}
N_VCM[n].xx=N_VCM[n].r11*N_VCM[n].crx+N_VCM[n].r12*N_VCM[n].cry+N_VCM[n].r13*N_VCM[n].crz;
N_VCM[n].yy=N_VCM[n].r21*N_VCM[n].crx+N_VCM[n].r22*N_VCM[n].cry+N_VCM[n].r23*N_VCM[n].crz;
N_VCM[n].zz=N_VCM[n].r31*N_VCM[n].crx+N_VCM[n].r32*N_VCM[n].cry+N_VCM[n].r33*N_VCM[n].crz;
}
if(OBSM==1)
{
for(i=0;i<Nsolw;i++)
{
pwl=ww+i;
n=pwl->pn;
if(n > 10000000)continue;
rbp=sqrt((ylc-pwl->y)*(ylc-pwl->y)+(zlc-pwl->z)*(zlc-pwl->z));
if(rbp >= Lr0)continue;
xl=pwl->xa-N_VCM[n].xcm;
yl=pwl->ya-N_VCM[n].ycm;
zl=pwl->za-N_VCM[n].zcm;
px=N_VCM[n].xx;
py=N_VCM[n].yy;
pz=N_VCM[n].zz;
vxs=crossx(px,py,pz,xl,yl,zl);
vys=crossy(px,py,pz,xl,yl,zl);
vzs=crossz(px,py,pz,xl,yl,zl);
pwl->vrx+=N_VCM[n].Vx-N_VCM[n].Vnx+vxs;
pwl->vry+=N_VCM[n].Vy-N_VCM[n].Vny+vys;
pwl->vrz+=N_VCM[n].Vz-N_VCM[n].Vnz+vzs;
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
xl=pwl->x;
yl=pwl->y;
zl=pwl->z;
rpc=sqrt((xl-rc->xcm)*(xl-rc->xcm)+(yl-rc->ycm)*(yl-rc->ycm)+(zl-rc->zcm)*(zl-rc->zcm));
if(rpc<Rrbc)
{
rc->Momx+= pwl->vrx-pwl->vx;
rc->Momy+= pwl->vry-pwl->vy;
rc->Momz+= pwl->vrz-pwl->vz;
rc->Amomx+=crossx((xl-rc->xcm),(yl-rc->ycm),(zl-rc->zcm),(pwl->vrx-pwl->vx),(pwl->vry-pwl->vy),(pwl->vrz-pwl->vz));
rc->Amomy+=crossy((xl-rc->xcm),(yl-rc->ycm),(zl-rc->zcm),(pwl->vrx-pwl->vx),(pwl->vry-pwl->vy),(pwl->vrz-pwl->vz));
rc->Amomz+=crossz((xl-rc->xcm),(yl-rc->ycm),(zl-rc->zcm),(pwl->vrx-pwl->vx),(pwl->vry-pwl->vy),(pwl->vrz-pwl->vz));
ic=Nrbcp+10;
}
}
}
ki=0;
kj=0;
hux[ki++] = 0.0;
hvx[kj++] = 0.0;
hax[0] = 0.0;
hbx[0] = 0.0;
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip != 1)continue;
if(rc->xcm <= xsl+Rrbc+0.55)
{
hux[ki++]=(double)rc->id;
hux[ki++]=rc->Momx;
hux[ki++]=rc->Momy;
hux[ki++]=rc->Momz;
hux[ki++]=rc->Amomx;
hux[ki++]=rc->Amomy;
hux[ki++]=rc->Amomz;
}
else if(rc->xcm >= xer-Rrbc-0.55)
{
hvx[kj++]=(double)rc->id;
hvx[kj++]=rc->Momx;
hvx[kj++]=rc->Momy;
hvx[kj++]=rc->Momz;
hvx[kj++]=rc->Amomx;
hvx[kj++]=rc->Amomy;
hvx[kj++]=rc->Amomz;
}
}
if(ki == 1) ki=0;
if(kj == 1) kj=0;
hux[0]=(double)ki;
hvx[0]=(double)kj;
if(numprocs > 1){
MPI_Sendrecv(&hux[0],ki,MPI_DOUBLE,left,195,&hax[0],5000,MPI_DOUBLE,right,195,comm1d,&status);
MPI_Sendrecv(&hvx[0],kj,MPI_DOUBLE,right,196,&hbx[0],5000,MPI_DOUBLE,left,196,comm1d,&status);
}
///*****************************************//
ki = (int)hax[0];
kj = 1;
while( kj < ki){
i=(int)hax[kj++];
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip != 1)continue;
if(rc->id == i)
{
rc->Momx+=hax[kj++];
rc->Momy+=hax[kj++];
rc->Momz+=hax[kj++];
rc->Amomx+=hax[kj++];
rc->Amomy+=hax[kj++];
rc->Amomz+=hax[kj++];
ic=Nrbcp+10;
}
}
if(ic < Nrbcp+3)
{
printf("collERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_ip");
printf(" rank=%i id=%i xsl=%i xer=%i \n",rank,i,xsl,xer);
exit(1);
}
}
ki = (int)hbx[0];
kj = 1;
while( kj < ki){
i=(int)hbx[kj++];
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip != 1)continue;
if(rc->id == i)
{
rc->Momx+=hbx[kj++];
rc->Momy+=hbx[kj++];
rc->Momz+=hbx[kj++];
rc->Amomx+=hbx[kj++];
rc->Amomy+=hbx[kj++];
rc->Amomz+=hbx[kj++];
ic=Nrbcp+10;
}
}
if(ic < Nrbcp+3)
{
printf("collERRRRRRRRRRROOOOOOOOOOOOOORRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRRR_ipbb");
printf(" rank=%i id=%i xsl=%i xer=%i \n",rank,i,xsl,xer);
exit(1);
}
}
for(ic=0;ic<Nrbcp;ic++)
{
rc=rbc+ic;
if(rc->ip == 2)continue;
rc->Vxcm+=(rc->Momx/Mrbc);
rc->Vycm+=(rc->Momy/Mrbc);
rc->Vzcm+=(rc->Momz/Mrbc);
rc->wx+=(rc->Amomx/Irbc);
rc->wy+=(rc->Amomy/Irbc);
rc->wz+=(rc->Amomz/Irbc);
rc->ex+=(crossx(rc->wx,rc->wy,rc->wz,rc->ex,rc->ey,rc->ez))*dt;
rc->ey+=(crossy(rc->wx,rc->wy,rc->wz,rc->ex,rc->ey,rc->ez))*dt;
rc->ez+=(crossz(rc->wx,rc->wy,rc->wz,rc->ex,rc->ey,rc->ez))*dt;
rc->ex=rc->ex/sqrt(rc->ex*rc->ex+rc->ey*rc->ey+rc->ez*rc->ez);
rc->ey=rc->ey/sqrt(rc->ex*rc->ex+rc->ey*rc->ey+rc->ez*rc->ez);
rc->ez=rc->ez/sqrt(rc->ex*rc->ex+rc->ey*rc->ey+rc->ez*rc->ez);
rc->Momx=0.0;
rc->Momy=0.0;
rc->Momz=0.0;
rc->Amomx=0.0;
rc->Amomy=0.0;
rc->Amomz=0.0;
}
}
ki=2;
hvx[0] = 0.0;
hvx[1] = 0.0;
hbx[0] = 0.0;
hbx[1] = 0.0;
for(i=0;i<Nsolpa;i++)
{
srd=ss+i;
n=srd->pn;
if(n > 10000000) continue;
xl=srd->xa-N_VCM[n].xcm;
yl=srd->ya-N_VCM[n].ycm;
zl=srd->za-N_VCM[n].zcm;
px=N_VCM[n].xx;
py=N_VCM[n].yy;
pz=N_VCM[n].zz;
vxs=crossx(px,py,pz,xl,yl,zl);
vys=crossy(px,py,pz,xl,yl,zl);
vzs=crossz(px,py,pz,xl,yl,zl);
srd->vx+=N_VCM[n].Vx-N_VCM[n].Vnx+vxs;
srd->vy+=N_VCM[n].Vy-N_VCM[n].Vny+vys;
srd->vz+=N_VCM[n].Vz-N_VCM[n].Vnz+vzs;
if(i >= Nsolp)
{
hvx[ki++]=(double)srd->n;
hvx[ki++]= srd->vx;
hvx[ki++]= srd->vy;
hvx[ki++]= srd->vz;
}
}
hvx[1]=(double)ki;
v2=0;
for(i=0;i<Naf;i++)
{
a=at+i;
n=a->pn;
if(n > 10000000) continue;
xl=a->xa-N_VCM[n].xcm;
yl=a->ya-N_VCM[n].ycm;
zl=a->za-N_VCM[n].zcm;
px=N_VCM[n].xx;
py=N_VCM[n].yy;
pz=N_VCM[n].zz;
vxs=crossx(px,py,pz,xl,yl,zl);
vys=crossy(px,py,pz,xl,yl,zl);
vzs=crossz(px,py,pz,xl,yl,zl);
a->vx+=N_VCM[n].Vx-N_VCM[n].Vnx+vxs;
a->vy+=N_VCM[n].Vy-N_VCM[n].Vny+vys;
a->vz+=N_VCM[n].Vz-N_VCM[n].Vnz+vzs;
v2+=a->vx*a->vx+a->vy*a->vy+a->vz*a->vz;
if(a->typ != 0)
{
if(a->xreal >= xsl && a->xreal < xer)
{
printf("\n rankkknnnnnn=%i n=%i typ=%i xl=%e xsl=%i xer=%i \n",rank,n,a->typ,a->xreal,xsl,xer);
printf("\n err_message(7)");
exit(1);
}
hvx[ki++]=(double)i;
hvx[ki++]= a->vx;
hvx[ki++]= a->vy;
hvx[ki++]= a->vz;
}
}
if(ki == 2) ki=0;
hvx[0]=(double)ki;
taf+=0.5*(tbf-v2*inv_Naf);
if(numprocs > 1){
if(shift1 > 0)
MPI_Sendrecv(&hvx[0],ki,MPI_DOUBLE,right,9,&hbx[0],600000,MPI_DOUBLE,left,9,comm1d,&status);
else
MPI_Sendrecv(&hvx[0],ki,MPI_DOUBLE,left,10,&hbx[0],600000,MPI_DOUBLE,right,10,comm1d,&status);
}
ki = (int)hbx[0];
j = (int)hbx[1];
kj = 2;
while( kj < ki){
if( kj < j){
i=(int)hbx[kj++];
srd=ss+i;
srd->vx=hbx[kj++];
srd->vy=hbx[kj++];
srd->vz=hbx[kj++];
}
else{
i=(int)hbx[kj++];
a=at+i;
a->vx=hbx[kj++];
a->vy=hbx[kj++];
a->vz=hbx[kj++];
}
}
}
void results()
{
unsigned long int i;
double xcm,ycm,zcm,r2,v2,pe,tfx,tfy,tfz,tfnor1,tfnor2,tfper,ttx,tty,ttz,tvx,tvy,tvz;
double tvnor1,tvnor2,tvper,tax,tay,taz,inv_d;
struct aftr *a;
struct sol *srd;
xcm=ycm=zcm=r2=v2=pe=tfx=tfy=tfz=tfnor1=tfnor2=tfper=tvnor1=tvnor2=tvper=ttx=tty=ttz=tvx=tvy=tvz=0;
tax=tay=taz=0;
for(i=0;i<Naf;i++)
{
a=at+i;
r2+=dist2(a->x0,a->xreal,a->y0,a->yreal,a->z0,a->zreal);
v2+=a->vx*a->vx+a->vy*a->vy+a->vz*a->vz;
pe+=a->PE;
tfx+=a->ax;
tfy+=a->ay;
tfz+=a->az;
tvx+=a->vx;
tvy+=a->vy;
tvz+=a->vz;
ttx+=crossx(a->xreal,a->yreal,a->zreal,a->ax,a->ay,a->az);
tty+=crossy(a->xreal,a->yreal,a->zreal,a->ax,a->ay,a->az);
ttz+=crossz(a->xreal,a->yreal,a->zreal,a->ax,a->ay,a->az);
inv_d=1./(a->xreal*a->xreal+a->yreal*a->yreal+a->zreal*a->zreal);
tax+=crossx(a->xreal,a->yreal,a->zreal,a->vx,a->vy,a->vz)*inv_d;
tay+=crossy(a->xreal,a->yreal,a->zreal,a->vx,a->vy,a->vz)*inv_d;
taz+=crossz(a->xreal,a->yreal,a->zreal,a->vx,a->vy,a->vz)*inv_d;
tfnor1+=(a->ax*a->nx+a->ay*a->ny+a->az*a->nz);
tfnor2+=(a->ax*a->tx+a->ay*a->ty+a->az*a->tz);
tfper+=(a->ax*(cnt+a->plane)->dirx+a->ay*(cnt+a->plane)->diry+a->az*(cnt+a->plane)->dirz);
tvnor1+=(a->vx*a->nx+a->vy*a->ny+a->vz*a->nz);
tvnor2+=(a->vx*a->tx+a->vy*a->ty+a->vz*a->tz);
tvper+=(a->vx*(cnt+a->plane)->dirx+a->vy*(cnt+a->plane)->diry+a->vz*(cnt+a->plane)->dirz);
a->PE=0;
}
//printf("Nsol=%le\n",v2/(double)(Nsol+Naf));
R2+=r2*inv_Naf;
R2CM+=dist2((cm+0)->xc0,(cm+0)->xreal,(cm+0)->yc0,(cm+0)->yreal,(cm+0)->zc0,(cm+0)->zreal);
KE+=0.5*v2*inv_Naf;
PE+=pe*inv_Naf;
totforx+=tfx*inv_Naf;
totfory+=tfy*inv_Naf;
totforz+=tfz*inv_Naf;
tottorx+=ttx*inv_Naf;
tottory+=tty*inv_Naf;
tottorz+=ttz*inv_Naf;
totangx+=tax*inv_Naf;
totangy+=tay*inv_Naf;
totangz+=taz*inv_Naf;
TNF1+=tfnor1*inv_Naf;
TNF2+=tfnor2*inv_Naf;
TPF+=tfper*inv_Naf;
TNV1+=tvnor1*inv_Naf;
TNV2+=tvnor2*inv_Naf;
TPV+=tvper*inv_Naf;
totvx+=tvx*inv_Naf;
totvy+=tvy*inv_Naf;
totvz+=tvz*inv_Naf;
//tvx=tvy=tvz=0;
for(i=0;i<Nsol;i++)
{
srd=ss+i;
v2+=srd->vx*srd->vx+srd->vy*srd->vy+srd->vz*srd->vz;
}
vsol2+=v2/(double)(Nsol+Naf);
}
void average()
{
double Nav,cp,cp1,e2edis,cp2,invcol;
Nav=1./(double)classwidth;
invcol=1./((double)classwidth*dtaf/dt);
//printf("sol=%le\n",vsol2/3.);
e2edis=sqrt(dist2((at+0)->xreal,(at+length-1)->xreal,(at+0)->yreal,(at+length-1)->yreal,(at+0)->zreal,(at+length-1)->zreal));
fp=fopen("diffusion.dat","a");
fprintf(fp,"%le\t%le\t%le\t",totdt,R2*Nav,R2CM*Nav);
fprintf(fp,"%le\t%le\t%le\t%le\t%le\n",KE*Nav,PE*Nav,(PE+KE)*Nav,vsol2*0.33333*Nav,e2edis*inv_tlength);
fclose(fp);
cp=(at+0)->nx*(at+0)->nx0+(at+0)->ny*(at+0)->ny0+(at+0)->nz*(at+0)->nz0;
cp1=(at+(Naf-1))->nx*(at+(Naf-1))->nx0+(at+(Naf-1))->ny*(at+(Naf-1))->ny0+(at+(Naf-1))->nz*(at+(Naf-1))->nz0;
cp2=(at+(Naf-1))->nx*(at+0)->nx+(at+(Naf-1))->ny*(at+0)->ny+(at+(Naf-1))->nz*(at+0)->nz;
fp=fopen("for-vel.dat","a");
fprintf(fp,"%le\t%le\t%le\t",totdt,totforx*Nav,totfory*Nav);
fprintf(fp,"%le\t%le\t%le\t%le\t",totforz*Nav,totvx*Nav,totvy*Nav,totvz*Nav);
fprintf(fp,"%le\t%le\t%le\t%le\t%le\t%le\n",tottorx*Nav,tottory*Nav,tottorz*Nav,totangx*Nav,totangy*Nav,totangz*Nav);
fclose(fp);
fp=fopen("angle.dat","a");
fprintf(fp,"%le\t%le\t%le\t%le\n",totdt,cp,cp1,cp2);
fclose(fp);
fp=fopen("fvpp.dat","a");
fprintf(fp,"%le\t%le\t%le\t%le\t%le\t%le\t%le\t%le\n",totdt,TNF1*Nav,TNF2*Nav,TPF*Nav,TNV1*Nav,TNV2*Nav,TPV*Nav,taf*invcol);
fclose(fp);
R2=R2CM=KE=PE=totforx=totfory=totforz=totvx=totvy=totvz=vsol2=TNF1=TNF2=TPF=TNV1=TNV2=TPV=tottorx=tottory=tottorz=taf=0;
totangx=totangy=totangz=0;
}
void transverse()
{
unsigned long int i,out,k,j,tl,out1;
double dam1x,dam1y,dam1z,d2,d,inv_d,fx,fy,fz,D,kh,strans;
double nx,ny,nz,ux,uy,uz,tx,ty,tz,ct,ct1,g;
struct aftr *ram2,*ram1;
i=out=0;
for(k=0;k<npoints;k++)
{
ram2=(at+(length*k));
j=out1=0;
tl=length-1;
do
{
if(j==0)
ram1=ram2->flagella[0];
else
ram1=ram2->flagella[1];/*a-1*/
nx=ram2->nx;
ny=ram2->ny;
nz=ram2->nz;
dam1x=(cnt+ram2->plane)->dirx;
dam1y=(cnt+ram2->plane)->diry;
dam1z=(cnt+ram2->plane)->dirz;
tx=crossx(nx,ny,nz,dam1x,dam1y,dam1z);
ty=crossy(nx,ny,nz,dam1x,dam1y,dam1z);
tz=crossz(nx,ny,nz,dam1x,dam1y,dam1z);
//inv_d=0;
d2=tx*tx+ty*ty+tz*tz;
d=sqrt(d2);
inv_d=1./d;
tx*=inv_d;
ty*=inv_d;
tz*=inv_d;
ram2->tx=tx;
ram2->ty=ty;
ram2->tz=tz;
//printf("trans1 %ld\t%ld\t%ld\n",ram2->n,ram1->n,j);
if(j==tl-1)
{
ram1->tx=tx;
ram1->ty=ty;
ram1->tz=tz;
//printf("trans11111 %ld\t%ld\t%ld\n",ram2->n,ram1->n,j);
}
//printf("bending body %ld\t%ld\t%ld\t%ld\t%ld\t\t%ld\t%le\n",ram2->n,ram1->n,ra->n,rap1->n,rap2->n,j,ra->cangle);
ram2=ram1;/*a-1*/
j++;
//ra->fz+=10.*sin(2.*PI*sqrt(caa));
if(ram1==(at+(k*length+tl)))/*(cm+0)->mtube[k])*/
out1=1;
}
while(out1==0);
}
}
void flalink(struct aftr *a,struct aftr *b)
{
unsigned long int i;
i=0;
while (a->flagella[i] && i<2)
i++;
if (i<2)
a->flagella[i]=b;
else
{
printf("MORE THAN 2 NEIGHB in flalink. !!! %ld\t%ld\n",a->flagella[0]->n,a->flagella[1]->n);
exit(1);
}
}
void hellink(struct aftr *a,struct aftr *b)
{
unsigned long int i;
i=0;
while (a->helix[i] && i<2)
i++;
if (i<2)
a->helix[i]=b;
else
{
printf("MORE THAN 2 NEIGHB in hellink. !!! %ld\t%ld\t%ld\n",a->n,a->helix[0]->n,a->helix[1]->n);
exit(1);
}
}
void velupdate(double nx,double ny,double nz,double tx1,double ty1,double tz1, unsigned long int afp, unsigned long int mn)
{
unsigned long int i,n;
double nr,tg,bi,vx,vy,vz,tx2,ty2,tz2,inv_d;
struct sol *srd;
srd=ss+mn;
nr=rayleighrandom();
tg=gaussianrandom();
bi=gaussianrandom();
tx2=crossx(nx,ny,nz,tx1,ty1,tz1);
ty2=crossy(nx,ny,nz,tx1,ty1,tz1);
tz2=crossz(nx,ny,nz,tx1,ty1,tz1);
inv_d=1./sqrt(tx2*tx2+ty2*ty2+tz2*tz2);
tx2*=inv_d;
ty2*=inv_d;
tz2*=inv_d;
vx=nx*nr+tx1*tg+tx2*bi;
vy=ny*nr+ty1*tg+ty2*bi;
vz=nz*nr+tz1*tg+tz2*bi;
vaddx+=(srd->vx-vx);
vaddy+=(srd->vy-vy);
vaddz+=(srd->vz-vz);
srd->vx=vx;
srd->vy=vy;
srd->vz=vz;
srd->in=1;
}
void velupdate2(unsigned long int afp, unsigned long int mn)
{
double vx,vy,vz;
struct aftr *a,*a1;
struct sol *srd;
srd=ss+mn;
a=at+afp;
vx=a->vx;
vy=a->vy;
vz=a->vz;
vaddx+=(srd->vx-vx);
vaddy+=(srd->vy-vy);
vaddz+=(srd->vz-vz);
/* a->vx+=(srd->vx-vx);
a->vy+=(srd->vy-vy);
a->vz+=(srd->vz-vz);
*/
srd->vx=vx;
srd->vy=vy;
srd->vz=vz;
srd->in=1;
}
double rayleighrandom()
{
return(SRTEMP*sqrt(-2.*log(1.-drand48())));
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment