c Computes moment of inertia tensor, volume, surface area, and integrated c curvatures (Gaussian and mean) include 'mpif.h' real rr(1000,1000),len,p(0:200,0:500) real xg(1000),wg(1000),phi,pi real theta,I11,I22,I33,I12,I13,I23 real I11s,I22s,I33s,I12s,I13s,I23s real volsave,volume,ggg,gggsave,curv,curvs real hh,kk,rr1,rr2,meanrc,theta1,phi1,theta2,phi2 real zzz,alphaL,betaL,xn,yn,zn real alphaW,betaW,wwwx,wwwy,wwwz real alphaT,betaT,tttx,ttty,tttz complex y(0:200,-250:250),a(0:200,-250:250) complex r,rtheta,rphi,rphiphi,rthth,rphith complex yphi(0:200,-250:250),sa,yphith(0:200,-250:250) complex ytheta(0:200,-250:250),f1,ythth(0:200,-250:250) complex yphiphi(0:200,-250:250),r1 complex ze,zf,zg,zl,zm,znn,zq,zh,zk,zq1 complex aaa,bbb,ccc,ddd,ddd2 complex thx,thy,thz,phix,phiy,phiz complex nxphi,nxth,nyphi,nyth,nzphi,nzth complex testhh,testkk,inttest integer pix(200,200,200),part(1000000,3),ijkmn integer nsave,nnn c needs to be changed with filename length changes character*26 filename,filename2 character*30 filename3 character*27 namef character*29 namef1 character*56 namef13 c do not need to be changed character*19 namef2 character*42 namef12 character*22 namef122 character*27 namef133 integer ierr,npes,me,Ntot,first,last,ios,first1 integer nnn1 integer n2,ind,ss,ss2 integer rem real length,width,thick,leth,with,thth real lengths,widths,thicks,leths,withs,thths real ssa,ssas,eqsa,eqsas,eqI,eqIs,trI,trIs,eqdiam,eqdiams c================================================================ c mpi initialization call mpi_init(ierr) c mpi structure : return number of processors and process number call mpi_comm_size(MPI_COMM_WORLD,npes, ierr) call mpi_comm_rank(MPI_COMM_WORLD,me,ierr) c================================================================== c set the number of anm files Ntot=14680 pi=4.0*atan(1.0) ng=120 c===================================================================== c create mpi vector data type call mpi_type_vector(1,ng,0, MPI_REAL8, ivectype, ierr) call mpi_type_commit(ivectype, ierr) c ============================================================== c use of mpi time to know how long each processor takes to run tc0 = MPI_Wtime() tc1 = MPI_Wtime() tar = tc1-tc0 tc0 = MPI_Wtime() c************************************************************* c ===== only the root processor reads the Gaussian file =========== if(me.eq.0) then open(unit=30,file='gauss120.dat') do 1100 i=1,ng read(30,*) xg(i),wg(i) 1100 continue close(30) endif c************************************************************* c root processor send a 2 vector data to all other processors c the vector data refers to those read from gaussien file call mpi_bcast (xg(1), 1, ivectype,0, MPI_COMM_WORLD, ierr) call mpi_bcast (wg(1), 1, ivectype, 0, MPI_COMM_WORLD, ierr) c ===================================================================================================== open(unit=29,file='anm.lis') c **** creation of Sample-V-a-* files which depend on the process rank***** Write (namef(20:23),'(I4)') 1000+me namef(1:20)='ASTMPM2-un-geom-len-' namef(24:27)='.dat' open(unit=12,file=namef) c Write (namef1(22:25),'(I4)') 1000+me c namef1(1:22)='ASTMPM2-un-SHrej-geom-' c namef1(26:29)='.dat' c open(unit=13,file=namef1) c calculation of the number of files that each processor must read last= Ntot/npes first1= Ntot/npes first = me*first1 + 1 rem = Ntot-npes*first1 if (rem.ne.0) then if(me.gt.(npes-rem)) then last=last+1 first=first+(me-npes-1+rem) end if end if c *** set the first file that the processor should read c Write (filename2(14:20),'(I7)') 1000000+first c filename2(1:14)='TypeI-gas-anm-' c filename2(21:24)='.dat' c if-test to find the right file in 'anm.lis' from where each processor should start 8888 continue read(29,2999) ijkmn,filename 2999 format(i6,1x,a26) c if(filename.ne.filename2) goto 8888 if(ijkmn.lt.first) goto 8888 c =========================================================================================== c the last proc has been set to read the rest of files if (me.eq.npes-1)then last=last+mod(Ntot,npes) endif c___________________________________________________________________________________________ do 3800 lmn=1,last filename3(1:4)='anm/' filename3(5:30)=filename open(unit=19,file=filename3) c open(unit=19,file=filename) c Read in expansion coefficients c nnn is how many y's were used in series (n limit) c read in stored coefficients for particle of interest c write(7,*) ' nnn = ',nnn do 390 n=0,26 do 390 m=n,-n,-1 read(19,*) ii,jj,a1,a2 a(n,m)=cmplx(a1,a2) 390 continue rewind(19) c initialize stored values in case no value of nnn will give c an acceptable value of the integrated Gaussian curvature nsave=0 gggsave=0.0 volsave=0.0 xmins=0.0 ymins=0.0 zmins=0.0 xmaxs=0.0 ymaxs=0.0 zmaxs=0.0 lengths=0.0 widths=0.0 thicks=0.0 leths=0.0 withs=0.0 thths=0.0 ssas=0.0 eqsas=1.0 eqdiams=0.0 curvs=0.0 trIs=0.0 eqIs=1.0 I11s=0.0 I22s=0.0 I33s=0.0 I23s=0.0 I13s=0.0 do 5555 npq=1,13 nnn=2*npq c Compute moment of inertia tensor and volume factor=pi*pi/2. I11=0.0 I22=0.0 I33=0.0 I12=0.0 I23=0.0 I13=0.0 volume=0.0 xmax=0.0 ymax=0.0 zmax=0.0 xmin=10000.0 ymin=10000.0 zmin=10000.0 do 430 i=1,ng do 430 j=1,ng theta=0.5*pi*xg(i)+0.5*pi phi=pi*xg(j)+pi call harm(theta,phi,nnn,p,y) r1=cmplx(0.0,0.0) do 460 n=0,nnn do 460 m=n,-n,-1 r1=r1+a(n,m)*y(n,m) 460 continue xx=real(r1)*sin(theta)*cos(phi) yy=real(r1)*sin(theta)*sin(phi) zz=real(r1)*cos(theta) if(xmax.lt.xx) xmax=xx if(ymax.lt.yy) ymax=yy if(zmax.lt.zz) zmax=zz if(xmin.gt.xx) xmin=xx if(ymin.gt.yy) ymin=yy if(zmin.gt.zz) zmin=zz I11=I11+0.2*wg(i)*wg(j)*(r1**5)*sin(theta)* & (1.-(cos(phi)*sin(theta))**2) I22=I22+0.2*wg(i)*wg(j)*(r1**5)*sin(theta)* & (1.-(sin(phi)*sin(theta))**2) I33=I33+0.2*wg(i)*wg(j)*(r1**5)*sin(theta)* & (1.-(cos(theta))**2) I12=I12-0.2*wg(i)*wg(j)*(r1**5)*sin(theta)* & cos(phi)*sin(phi)*(sin(theta)**2) I13=I13-0.2*wg(i)*wg(j)*(r1**5)*sin(theta)* & cos(phi)*cos(theta)*sin(theta) I23=I23-0.2*wg(i)*wg(j)*(r1**5)*sin(theta)* & sin(phi)*cos(theta)*sin(theta) volume=volume+wg(i)*wg(j)*(r1**3)*sin(theta)/3. 430 continue volume=volume*factor I11=I11*factor I22=I22*factor I33=I33*factor I12=I12*factor I13=I13*factor I23=I23*factor I11=I11/volume I22=I22/volume I33=I33/volume I12=I12/volume I13=I13/volume I23=I23/volume c call flush(7) c compute surface area and integrated curvatures c using spherical harmonic expansion factor=pi*pi/2. sa=cmplx(0.,0.) zh=cmplx(0.,0.) zk=cmplx(0.,0.) meanrc=0. inttest=cmplx(0.,0.) do 420 i=1,ng do 420 j=1,ng theta=0.5*pi*xg(i)+0.5*pi phi=pi*xg(j)+pi call harm(theta,phi,nnn+2,p,y) call surharm(theta,phi,nnn,p,y,yphi,ytheta,yphiphi,ythth, & yphith) r=a(0,0)*y(0,0) rtheta=cmplx(0.,0.) rthth=cmplx(0.,0.) rphi=cmplx(0.,0.) rphiphi=cmplx(0.,0.) rphith=cmplx(0.,0.) c derivatives of (0,0) term is zero, since it is a constant do 22 n=1,nnn do 22 m=n,-n,-1 r=r+a(n,m)*y(n,m) rphi=rphi+a(n,m)*yphi(n,m) rphiphi=rphiphi+a(n,m)*yphiphi(n,m) rphith=rphith+a(n,m)*yphith(n,m) rtheta=rtheta+a(n,m)*ytheta(n,m) rthth=rthth+a(n,m)*ythth(n,m) 22 continue ct=cos(theta) st=sin(theta) cp=cos(phi) sp=sin(phi) ddd=rphi**2+rtheta*rtheta*st*st+r*r*st*st ddd2=csqrt(ddd) f1=r*ddd2 ze=rtheta**2+r**2 zf=rtheta*rphi zg=rphi**2+r*r*st*st c derivatives of coordinate thx=(rtheta*st+r*ct)*cp thy=(rtheta*st+r*ct)*sp thz=rtheta*ct-r*st phix=st*(rphi*cp-r*sp) phiy=st*(rphi*sp+r*cp) phiz=rphi*ct c normal(x), phi derivative aaa=(rphi**2)*sp+r*rphiphi*sp+r*rphi*cp & -rphi*rtheta*st*ct*cp-r*rphith*st*ct*cp & +r*rtheta*st*ct*sp+2.*r*rphi*st*st*cp & -r*r*st*st*sp bbb=r*rphi*sp-r*rtheta*st*ct*cp+r*r*st*st*cp ccc=rphi/r+(rphi*rphiphi+rtheta*rphith*st*st+r*rphi*st*st) & /ddd nxphi=(1./r/ddd2)*(aaa-bbb*ccc) c normal(y), phi derivative aaa=-rphi*rphi*cp-r*rphiphi*cp+r*rphi*sp-rphi*rtheta*st*ct*sp & -r*rphith*st*ct*sp-r*rtheta*st*ct*cp+2.*r*rphi*st*st*sp & +r*r*st*st*cp bbb=-r*rphi*cp-r*rtheta*st*ct*sp+r*r*st*st*sp nyphi=(1./r/ddd2)*(aaa-bbb*ccc) c normal(z), phi derivative aaa=rphi*rtheta*st*st+r*rphith*st*st+2.*r*rphi*ct*st bbb=r*rtheta*st*st+r*r*ct*st nzphi=(1./r/ddd2)*(aaa-bbb*ccc) c normal(x), theta derivative aaa=rtheta*rphi*sp+r*rphith*sp-rtheta*rtheta*st*ct*cp & -r*rthth*st*ct*cp-r*rtheta*ct*ct*cp+r*rtheta*st*st*cp & +2.*r*rtheta*st*st*cp+2.*r*r*st*ct*cp bbb=r*rphi*sp-r*rtheta*st*ct*cp+r*r*st*st*cp ccc=(rphi*rphith+rtheta*rthth*st*st+r*rtheta*st*st & +rtheta*rtheta*st*ct+r*r*st*ct)/ddd+rtheta/r nxth=(1./r/ddd2)*(aaa-bbb*ccc) c normal(y), theta derivative aaa=-rtheta*rphi*cp-r*rphith*cp-rtheta*rtheta*st*ct*sp & -r*rthth*st*ct*sp-r*rtheta*ct*ct*sp+r*rtheta*st*st*sp & +2.*r*rtheta*st*st*sp+2.*r*r*st*ct*sp bbb=-r*rphi*cp-r*rtheta*st*ct*sp+r*r*st*st*sp nyth=(1./r/ddd2)*(aaa-bbb*ccc) c normal(z), theta derivative aaa=rtheta*rtheta*st*st+r*rthth*st*st+2.*r*rtheta*st*ct & +2.*r*rtheta*ct*st+r*r*ct*ct-r*r*st*st bbb=r*rtheta*st*st+r*r*ct*st nzth=(1./r/ddd2)*(aaa-bbb*ccc) zl=-(thx*nxth+thy*nyth+thz*nzth) znn=-(phix*nxphi+phiy*nyphi+phiz*nzphi) zq=-(phix*nxth+phiy*nyth+phiz*nzth) zq1=-(thx*nxphi+thy*nyphi+thz*nzphi) zm=-0.5*(zq1+zq) c integrated quantities inttest=inttest+wg(i)*wg(j)*f1 sa=sa+wg(i)*wg(j)*f1 zh=zh+wg(i)*wg(j)*f1*(ze*znn+zg*zl-2.*zf*zm)/2./(ze*zg-zf*zf) zk=zk+wg(i)*wg(j)*f1*(zl*znn-zm*zm)/(ze*zg-zf*zf) testhh=(ze*znn+zg*zl-2.*zf*zm)/2./(ze*zg-zf*zf) testkk=(zl*znn-zm*zm)/(ze*zg-zf*zf) hh=real((ze*znn+zg*zl-2.*zf*zm)/2./(ze*zg-zf*zf)) kk=real((zl*znn-zm*zm)/(ze*zg-zf*zf)) xk1=(hh+sqrt(hh*hh-kk)) xk2=(hh-sqrt(hh*hh-kk)) rr1=1./xk1 rr2=1./xk2 rr1=-0.5*(rr1+rr2) meanrc=meanrc+wg(i)*wg(j)*real(f1)*rr1 420 continue inttest=inttest*factor sa=sa*factor zh=zh*factor/sa zk=zk*factor meanrc=meanrc*factor/sa inttest=inttest/sa ggg=real(zk/4./pi) ssa=real(sa) rc=0.5*(6.*volume/pi)**(1./3.) eqsa=4.*pi*rc*rc eqdiam=2.*rc curv=-1./real(zh) trI=(I11+I22+I33)/3. eqI=0.4*rc*rc c lwt subroutine c call lwt(filename,nnn,a,me,length,width,thick,leth,with,thth) call lwt(filename,nnn,a,me,length,width,thick,leth,with & ,thth,xn,yn,zn,alphaL,betaL,wwwx,wwwy,wwwz,alphaW,betaW, & tttx,ttty,tttz,alphaT,betaT) c at this point, all geometrical parameters have been computed if(abs(ggg-1.0).lt.0.05) then c save all values if nnn > 10 if(nnn.gt.10)then nsave=nnn gggsave=ggg volsave=volume xmins=xmin ymins=ymin zmins=zmin xmaxs=xmax ymaxs=ymax zmaxs=zmax lengths=length widths=width thicks=thick leths=leth withs=with thths=thth ssas=ssa eqsas=eqsa eqdiams=eqdiam curvs=curv trIs=trI eqIs=eqI I11s=I11 I22s=I22 I33s=I33 I12s=I12 I23s=I23 I13s=I13 end if end if 5555 continue close(19) c if obtained abs(ggg-1) < 0.05 at nsave > 10, then it is a good c particle and write out data in geom file c use saved data - this will have been saved at maximum value of nsave c becuase of the do 3800 loop, nsave will be a its maximum value, c I don't want to distinguish between values of ggg, just as long they c are within 5 % of unity if(nsave.gt.10) then write(12,399) filename,xmins,xmaxs,ymins,ymaxs,zmins,zmaxs,volsave & ,ssas,ssas/eqsas,eqdiams,trIs/eqIs,nsave,gggsave,curvs,lengths & ,widths,thicks,leths,withs,thths,I11s,I22s,I33s,I12s,I13s,I23s & ,xn,yn,zn,alphaL,betaL,wwwx,wwwy,wwwz,alphaW,betaW & ,tttx,ttty,tttz,alphaT,betaT 399 format(a26,6f13.4,f18.6,f18.6,f8.3,f13.6,f9.3,i6,f10.5,f12.6 & ,6f14.6,6f15.7,15f13.8) call flush(12) end if c if nsave < 10 (could be zero if gg was never close to 1), then it is a bad c particle and write out data in SHrej file c use saved data - this will have been saved at maximum value of nsave c if(nsave.le.10) then c write(13,399) filename,xmins,xmaxs,ymins,ymaxs,zmins,zmaxs,volsave c & ,ssas,ssas/eqsas,eqdiams,trIs/eqIs,nsave,gggsave,curvs,lengths c & ,widths,thicks,leths,withs,thths,I11s,I22s,I33s,I12s,I13s,I23s c & ,xn,yn,zn,alphaL,betaL,wwwx,wwwy,wwwz,alphaW,betaW c & ,tttx,ttty,tttz,alphaT,betaT c end if c read the next filename in 'anm.lis' read(29,2999,iostat=ios) ijkmn,filename if(ios.ne.0) goto 1001 !if-test to exit at the end of the 'anm.lis' file 3800 continue c ____________________________________________________________________________________________ 1001 close(12) close(29) c calculation of time that each processor has done tc1 = MPI_Wtime() sec = (tc1-tc0) - tar print*,'proc = ',me,' finished on = ',sec,' seconds ' c mpi barrier used here to make sure that all processor should finish to continue call mpi_barrier ( MPI_COMM_WORLD, ierr) c ======================================================== c concatenate the results files if(me.eq.0)then do 2201 i=0,npes-1 Write (namef13(24:27),'(I4)') 1000+i namef13(1:24)='cat ASTMPM2-un-geom-len-' namef13(28:56)='.dat>>ASTMPM2-un-geom-len.dat' c Write (namef12(16:19),'(I4)') 1000+i c namef12(1:16)='cat length-data-' c namef12(20:42)='.dat >> length-data.dat' call system(namef13) c call system(namef12) 2201 continue c delete the temporary results files do 1101 i=0,npes-1 Write (namef122(15:18),'(I4)') 1000+i namef122(1:15)='rm length-data-' namef122(19:22)='.dat' c Write (namef133(20:23),'(I4)') 1000+i c namef133(1:20)='rm ArdalA-3-29-2013-geom-len-' c namef133(24:27)='.dat' call system(namef122) c call system(namef133) 1101 continue endif c=================================================== call mpi_finalize(ierr) end c ====================================================================================================== subroutine lwt(filename,nnn,a,me,length,width,thick,leth,with & ,thth,xn,yn,zn,alphaL,betaL,wwwx,wwwy,wwwz,alphaW,betaW, & tttx,ttty,tttz,alphaT,betaT) c Computes length (maximum distance between two points on rock), then maximum c distance between two points that is perpendicular to length, and then thickness, c which is the maximum distance perpendicular to both length and thickness. real rr(1000,1000),len,p(0:200,0:500) real xg(1000),wg(1000),phi,pi real theta,I11,I22,I33,I12,I13,I23 real hh,kk,rr1,rr2,meanrc real zzz,alphaL,betaL,xn,yn,zn real alphaW,betaW,wwwx,wwwy,wwwz real alphaT,betaT,tttx,ttty,tttz real, intent(out)::length,width,thick,leth,with,thth complex y(0:200,-250:250),a(0:200,-250:250) complex y1(0:200,-250:250),y2(0:200,-250:250) complex r,rtheta,rphi,rphiphi,rthth,rphith complex ze,zf,zg,zl,zm,zq,zh,zk,zq1 complex aaa,bbb,ccc,ddd,ddd2,r1,r2 complex thx,thy,thz,phix,phiy,phiz complex nxphi,nxth,nyphi,nyth,nzphi,nzth complex testhh,testkk,inttest integer pix(200,200,200),part(1000000,3) character*36 filename character*22 namef1 character*19 namef2 integer me c ******************************** c **** create the name of results files which depend on the procesor's rank c Write (namef1(15:18),'(I4)') 1000+me c namef1(1:15)='length-Sample-V-a-' c namef1(19:22)='.dat' Write (namef2(12:15),'(I4)') 1000+me namef2(1:12)='length-data-' namef2(16:19)='.dat' c ******************************** c open the files and start writing at the end of file thanks to the position option c it starts writing at the end to not erase what has written before ( at the first call of the routine) c open(unit=16,file=namef1,position='append') open(unit=15,file=namef2,position='append') pi=4.0*atan(1.0) c Read in expansion coefficients c nnn is how many y's were used in series (n limit) c read in stored coefficients for particle of interest nang=20 rmax=0.0 do 5555 i1=1,nang do 5555 j1=1,nang do 5555 i2=1,nang do 5555 j2=1,nang c Compute length theta1=0.001d0+(i1-1)*pi/float(nang) phi1=(j1-1)*2.*pi/float(nang) theta2=0.001d0+(i2-1)*pi/float(nang) phi2=(j2-1)*2.*pi/float(nang) call harm(theta1,phi1,nnn,p,y1) call harm(theta2,phi2,nnn,p,y2) r1=cmplx(0.0,0.0) r2=cmplx(0.0,0.0) do 460 n=0,nnn do 460 m=n,-n,-1 r1=r1+a(n,m)*y1(n,m) r2=r2+a(n,m)*y2(n,m) 460 continue xx1=real(r1)*sin(theta1)*cos(phi1) yy1=real(r1)*sin(theta1)*sin(phi1) zz1=real(r1)*cos(theta1) xx2=real(r2)*sin(theta2)*cos(phi2) yy2=real(r2)*sin(theta2)*sin(phi2) zz2=real(r2)*cos(theta2) rm=sqrt((xx1-xx2)**2+(yy1-yy2)**2+(zz1-zz2)**2) if(rmax.lt.rm) then rmax=rm t1=theta1 pp1=phi1 t2=theta2 pp2=phi2 end if 5555 continue length=rmax write(15,299) filename 299 format(a26) write(15,*) 'length = ',length write(15,*) 'theta1= ',t1,' phi1 = ',pp1 write(15,*) 'theta2= ',t2,' phi2 = ',pp2 write(15,*) call flush(15) do 5556 i1=1,nang-1 do 5556 j1=1,nang-1 do 5556 i2=1,nang-1 do 5556 j2=1,nang-1 c Compute refined length theta1=t1+(i1-10)*pi/float(nang)/10. phi1=pp1+(j1-10)*2.*pi/float(nang)/10. theta2=t2+(i2-10)*pi/float(nang)/10. phi2=pp2+(j2-10)*2.*pi/float(nang)/10. call harm(theta1,phi1,nnn,p,y1) call harm(theta2,phi2,nnn,p,y2) r1=cmplx(0.0,0.0) r2=cmplx(0.0,0.0) do 461 n=0,nnn do 461 m=n,-n,-1 r1=r1+a(n,m)*y1(n,m) r2=r2+a(n,m)*y2(n,m) 461 continue xx1=real(r1)*sin(theta1)*cos(phi1) yy1=real(r1)*sin(theta1)*sin(phi1) zz1=real(r1)*cos(theta1) xx2=real(r2)*sin(theta2)*cos(phi2) yy2=real(r2)*sin(theta2)*sin(phi2) zz2=real(r2)*cos(theta2) rm=sqrt((xx1-xx2)**2+(yy1-yy2)**2+(zz1-zz2)**2) if(rmax.lt.rm) then rmax=rm tt1=theta1 ppp1=phi1 tt2=theta2 ppp2=phi2 c compute unit vector along length xn=(xx2-xx1)/rmax yn=(yy2-yy1)/rmax zn=(zz2-zz1)/rmax end if 5556 continue length=rmax write(15,*) 'refined length = ',length write(15,*) 'theta1= ',tt1,' phi1 = ',ppp1 write(15,*) 'theta2= ',tt2,' phi2 = ',ppp2 write(15,*) call flush(15) c (xn,yn,zn) is unit vector along length rmax=0.0 c now compute width do 5557 i1=1,nang do 5557 j1=1,nang do 5557 i2=1,nang do 5557 j2=1,nang theta1=0.001+(i1-1)*pi/float(nang) phi1=(j1-1)*2.*pi/float(nang) theta2=0.001+(i2-1)*pi/float(nang) phi2=(j2-1)*2.*pi/float(nang) call harm(theta1,phi1,nnn,p,y1) call harm(theta2,phi2,nnn,p,y2) r1=cmplx(0.0,0.0) r2=cmplx(0.0,0.0) do 462 n=0,nnn do 462 m=n,-n,-1 r1=r1+a(n,m)*y1(n,m) r2=r2+a(n,m)*y2(n,m) 462 continue xx1=real(r1)*sin(theta1)*cos(phi1) yy1=real(r1)*sin(theta1)*sin(phi1) zz1=real(r1)*cos(theta1) xx2=real(r2)*sin(theta2)*cos(phi2) yy2=real(r2)*sin(theta2)*sin(phi2) zz2=real(r2)*cos(theta2) rm=sqrt((xx1-xx2)**2+(yy1-yy2)**2+(zz1-zz2)**2) c compute unit vector along trial width wx=(xx2-xx1)/rm wy=(yy2-yy1)/rm wz=(zz2-zz1)/rm dot=wx*xn+wy*yn+wz*zn if(abs(dot).le.0.1) then if(rmax.lt.rm) then rmax=rm t1=theta1 pp1=phi1 t2=theta2 pp2=phi2 wwx=wx wwy=wy wwz=wz end if end if 5557 continue width=rmax write(15,*) 'width = ',width write(15,*) 'theta1= ',t1,' phi1 = ',pp1 write(15,*) 'theta2= ',t2,' phi2 = ',pp2 write(15,*) ' length unit vector ',xn,yn,zn write(15,*) ' width unit vector ',wwx,wwy,wwz ddot=wwx*xn+wwy*yn+wwz*zn dot=acos(ddot) write(15,*) ' angle between lw = ',dot*180./pi-90.,' deg' write(15,*) call flush(15) c set final unit vector and angles for width to first iteration of width wwwx=wwx wwwy=wwy wwwz=wwz wtt1=t1 wppp1=pp1 wtt2=t2 wppp2=pp2 c refine width do 5558 i1=1,nang-1 do 5558 j1=1,nang-1 do 5558 i2=1,nang-1 do 5558 j2=1,nang-1 c Compute refined width theta1=t1+(i1-10)*pi/float(nang)/10. phi1=pp1+(j1-10)*2.*pi/float(nang)/10. theta2=t2+(i2-10)*pi/float(nang)/10. phi2=pp2+(j2-10)*2.*pi/float(nang)/10. call harm(theta1,phi1,nnn,p,y1) call harm(theta2,phi2,nnn,p,y2) r1=cmplx(0.0,0.0) r2=cmplx(0.0,0.0) do 463 n=0,nnn do 463 m=n,-n,-1 r1=r1+a(n,m)*y1(n,m) r2=r2+a(n,m)*y2(n,m) 463 continue xx1=real(r1)*sin(theta1)*cos(phi1) yy1=real(r1)*sin(theta1)*sin(phi1) zz1=real(r1)*cos(theta1) xx2=real(r2)*sin(theta2)*cos(phi2) yy2=real(r2)*sin(theta2)*sin(phi2) zz2=real(r2)*cos(theta2) rm=sqrt((xx1-xx2)**2+(yy1-yy2)**2+(zz1-zz2)**2) c compute unit vector along trial width wx=(xx2-xx1)/rm wy=(yy2-yy1)/rm wz=(zz2-zz1)/rm dot=wx*xn+wy*yn+wz*zn if(abs(dot).le.abs(ddot)) then if(rmax.lt.rm) then rmax=rm wtt1=theta1 wppp1=phi1 wtt2=theta2 wppp2=phi2 wwwx=wx wwwy=wy wwwz=wz end if end if 5558 continue width=rmax write(15,*) 'refined width = ',width write(15,*) 'theta1= ',wtt1,' phi1 = ',wppp1 write(15,*) 'theta2= ',wtt2,' phi2 = ',wppp2 write(15,*) ' length unit vector ',xn,yn,zn write(15,*) ' width unit vector ',wwwx,wwwy,wwwz dot=wwwx*xn+wwwy*yn+wwwz*zn dot=acos(dot) write(15,*) ' angle between lw = ',dot*180./pi-90.,' deg' write(15,*) call flush(15) rmax=0.0 c thickness: (xn,yn,zn) and (wwwx,wwwy,wwwz) are length and width unit vectors do 5559 i1=1,nang do 5559 j1=1,nang do 5559 i2=1,nang do 5559 j2=1,nang c Compute thickness theta1=0.001+(i1-1)*pi/float(nang) phi1=(j1-1)*2.*pi/float(nang) theta2=0.001+(i2-1)*pi/float(nang) phi2=(j2-1)*2.*pi/float(nang) call harm(theta1,phi1,nnn,p,y1) call harm(theta2,phi2,nnn,p,y2) r1=cmplx(0.0,0.0) r2=cmplx(0.0,0.0) do 464 n=0,nnn do 464 m=n,-n,-1 r1=r1+a(n,m)*y1(n,m) r2=r2+a(n,m)*y2(n,m) 464 continue xx1=real(r1)*sin(theta1)*cos(phi1) yy1=real(r1)*sin(theta1)*sin(phi1) zz1=real(r1)*cos(theta1) xx2=real(r2)*sin(theta2)*cos(phi2) yy2=real(r2)*sin(theta2)*sin(phi2) zz2=real(r2)*cos(theta2) rm=sqrt((xx1-xx2)**2+(yy1-yy2)**2+(zz1-zz2)**2) c compute unit vector along trial thickness tx=(xx2-xx1)/rm ty=(yy2-yy1)/rm tz=(zz2-zz1)/rm dotl=tx*xn+ty*yn+tz*zn dotw=wwwx*tx+wwwy*ty+wwwz*tz if(abs(dotl).le.0.1.and.abs(dotw).le.0.1) then if(rmax.lt.rm) then rmax=rm t1t=theta1 pp1t=phi1 t2t=theta2 pp2t=phi2 ttx=tx tty=ty ttz=tz end if end if 5559 continue thick=rmax write(15,*) 'thickness = ',thick write(15,*) 'theta1= ',t1t,' phi1 = ',pp1t write(15,*) 'theta2= ',t2t,' phi2 = ',pp2t write(15,*) ' length unit vector ',xn,yn,zn write(15,*) ' width unit vector ',wwwx,wwwy,wwwz write(15,*) ' thickness unit vector ',ttx,tty,ttz dotl=ttx*xn+tty*yn+ttz*zn dotw=wwwx*ttx+wwwy*tty+wwwz*ttz ddotl=abs(ttx*xn+tty*yn+ttz*zn) ddotw=abs(wwwx*ttx+wwwy*tty+wwwz*ttz) dotl=acos(dotl) dotw=acos(dotw) write(15,*) ' angle between lt = ',dotl*180./pi-90.,' deg' write(15,*) ' angle between wt = ',dotw*180./pi-90.,' deg' write(15,*) call flush(15) c set final unit vector and angles for thickness to c first iteration of thickness tttx=ttx ttty=tty tttz=ttz tt1t=t1t ppp1t=pp1t tt2t=t2t ppp2t=pp2t c refine thickness do 5560 i1=1,nang-1 do 5560 j1=1,nang-1 do 5560 i2=1,nang-1 do 5560 j2=1,nang-1 c Compute refined width theta1=t1t+(i1-10)*pi/float(nang)/10. phi1=pp1t+(j1-10)*2.*pi/float(nang)/10. theta2=t2t+(i2-10)*pi/float(nang)/10. phi2=pp2t+(j2-10)*2.*pi/float(nang)/10. call harm(theta1,phi1,nnn,p,y1) call harm(theta2,phi2,nnn,p,y2) r1=cmplx(0.0,0.0) r2=cmplx(0.0,0.0) do 465 n=0,nnn do 465 m=n,-n,-1 r1=r1+a(n,m)*y1(n,m) r2=r2+a(n,m)*y2(n,m) 465 continue xx1=real(r1)*sin(theta1)*cos(phi1) yy1=real(r1)*sin(theta1)*sin(phi1) zz1=real(r1)*cos(theta1) xx2=real(r2)*sin(theta2)*cos(phi2) yy2=real(r2)*sin(theta2)*sin(phi2) zz2=real(r2)*cos(theta2) rm=sqrt((xx1-xx2)**2+(yy1-yy2)**2+(zz1-zz2)**2) c compute unit vector along trial thickness tx=(xx2-xx1)/rm ty=(yy2-yy1)/rm tz=(zz2-zz1)/rm dotl=tx*xn+ty*yn+tz*zn dotw=wwwx*tx+wwwy*ty+wwwz*tz c both angles must decrease in order to even check if distance increases if(abs(dotl).le.ddotl.and.abs(dotw).le.ddotw) then if(rmax.lt.rm) then rmax=rm tt1t=theta1 ppp1t=phi1 tt2t=theta2 ppp2t=phi2 tttx=tx ttty=ty tttz=tz end if end if 5560 continue thick=rmax write(15,*) 'refined thickness = ',thick write(15,*) 'theta1= ',tt1t,' phi1 = ',ppp1t write(15,*) 'theta2= ',tt2t,' phi2 = ',ppp2t write(15,*) ' length unit vector ',xn,yn,zn write(15,*) ' width unit vector ',wwwx,wwwy,wwwz write(15,*) ' thickness unit vector ',tttx,ttty,tttz dotl=tttx*xn+ttty*yn+tttz*zn dotw=wwwx*tttx+wwwy*ttty+wwwz*tttz dotl=acos(dotl) dotw=acos(dotw) write(15,*) ' angle between lt = ',dotl*180./pi-90.,' deg' write(15,*) ' angle between wt = ',dotw*180./pi-90.,' deg' write(15,*) call flush(15) leth = length/thick with= width/thick thth= thick/thick c write(16,333) filename,length,width,thick, c & length/thick,width/thick,thick/thick c333 format(a24,6f11.6) c call flush(16) close(15) c close(16) c now work out Euler angles for each of these unit vectors (alpha,beta) c L zzz=sqrt(xn**2+yn**2) alphaL=acos(xn/zzz) if(yn.lt.0.0) alphaL=2.*pi-alphaL betaL=acos(zn) c W zzz=sqrt(wwwx**2+wwwy**2) alphaW=acos(wwwx/zzz) if(wwwy.lt.0.0) alphaW=2.*pi-alphaW betaW=acos(wwwz) c T zzz=sqrt(tttx**2+ttty**2) alphaT=acos(tttx/zzz) if(ttty.lt.0.0) alphaT=2.*pi-alphaT betaT=acos(tttz) return end c this is the factorial function function fac(j) integer j if(j.eq.0) then fac=1. return end if if(j.eq.1) then fac=1. return end if fac=1. do 10 i=1,j fac=fac*float(i) 10 continue return end subroutine harm(theta,phi,nnn,p,y) c computes spherical harmonics (complex) c for a value of x = cos(theta), phi=angle phi c so -1 < x < 1, P(n,m), -n < m < n, 0 < n . c uses two recursion relations, plus exact formulas c for the associated Legendre functions up to n=8 real p(0:200,0:500),x,s,fac complex y(0:200,-250:250) real phi,theta,pi,xn,xm pi=4.0*atan(1.0) x=cos(theta) s=sqrt(1.-x*x) do 50 n=0,nnn do 50 m=0,n p(n,m)=0.0 50 continue p(0,0)=1.0 p(1,0)=x p(1,1)=s p(2,0)=0.5*(3.*x*x-1.) p(2,1)=3.*x*s p(2,2)=3.*(1.-x*x) p(3,0)=0.5*x*(5.*x*x-3.) p(3,1)=1.5*(5.*x*x-1.)*s p(3,2)=15.*x*(1.-x*x) p(3,3)=15.*(s**3) p(4,0)=0.125*(35.*(x**4)-30.*x*x+3.) p(4,1)=2.5*(7.*x*x*x-3.*x)*s p(4,2)=7.5*(7.*x*x-1.)*(1.-x*x) p(4,3)=105.*x*(s**3) p(4,4)=105.*((1.-x*x)**2) p(5,0)=0.125*x*(63.*(x**4)-70.*x*x+15.) p(5,1)=0.125*15.*s*(21.*(x**4)-14.*x*x+1.) p(5,2)=0.5*105.*x*(1.-x*x)*(3.*x*x-1.) p(5,3)=0.5*105.*(s**3)*(9.*x*x-1.) p(5,4)=945.*x*((1.-x*x)**2) p(5,5)=945.*(s**5) p(6,0)=0.0625*(231.*(x**6)-315.*(x**4)+105.*x*x-5.) p(6,1)=0.125*21.*x*(33.*(x**4)-30.*x*x+5.)*s p(6,2)=0.125*105.*(1.-x*x)*(33.*(x**4)-18.*x*x+1.) p(6,3)=0.5*315.*(11.*x*x-3.)*x*(s**3) p(6,4)=0.5*945.*(1.-x*x)*(1.-x*x)*(11.*x*x-1.) p(6,5)=10395.*x*(s**5) p(6,6)=10395.*((1.-x*x)**3) p(7,0)=0.0625*x*(429.*(x**6)-693.*(x**4)+315.*x*x-35.) p(7,1)=0.0625*7.*s*(429.*(x**6)-495.*(x**4)+135.*x*x-5.) p(7,2)=0.125*63.*x*(1.-x*x)*(143.*(x**4)-110.*x*x+15.) p(7,3)=0.125*315.*(s**3)*(143.*(x**4)-66.*x*x+3.) p(7,4)=0.5*3465.*x*(1.-x*x)*(1.-x*x)*(13.*x*x-3.) p(7,5)=0.5*10395.*(s**5)*(13.*x*x-1.) p(7,6)=135135.*x*(1.-x*x)*(1.-x*x)*(1.-x*x) p(7,7)=135135.*(s**7) p(8,0)=(1./128.)*(6435.*(x**8)-12012.*(x**6)+6930.*(x**4)- & 1260.*x*x+35.) p(8,1)=0.0625*9.*x*s*(715.*(x**6)-1001.*(x**4)+385.*x*x-35.) p(8,2)=0.0625*315.*(1.-x*x)*(143.*(x**6)-143.*(x**4)+33.*x*x-1.) p(8,3)=0.125*3465.*x*(s**3)*(39.*(x**4)-26.*x*x+3.) p(8,4)=0.125*10395.*(1.-x*x)*(1.-x*x)*(65.*(x**4)-26.*x*x+1.) p(8,5)=0.5*135135.*x*(s**5)*(5.*x*x-1.) p(8,6)=0.5*135135.*((1.-x*x)**3)*(15.*x*x-1.) p(8,7)=2027025.*x*(s**7) p(8,8)=2027025.*((1.-x*x)**4) c now generate spherical harmonics for n=0,8 (follows Arfken) do 300 n=0,8 c does n=0 separately if(n.eq.0) then y(0,0)=1./sqrt(4.*pi) goto 300 end if do 200 m=n,-n,-1 if(m.ge.0) then y(n,m)=((-1.)**m)*sqrt( ((2*n+1)/4./pi)*fac(n-m)/fac(n+m) ) & *p(n,m)*cmplx(cos(m*phi),sin(m*phi)) goto 200 end if if(m.lt.0) then mm=-m y(n,m)=((-1.)**mm)*cmplx(real(y(n,mm)),-aimag(y(n,mm))) goto 200 end if 200 continue 300 continue c use recursion relations for n >= 9 c Do recursion on spherical harmonics, as they are better behaved numerically do 100 n=9,nnn do 99 m=0,n-2 xn=float(n-1) xm=float(m) y(n,m)=(2.*xn+1.)*x*y(n-1,m)-sqrt( (2.*xn+1.)*(xn**2-xm**2) & /(2.*xn-1.) )*y(n-2,m) y(n,m)=y(n,m)/sqrt( (2.*xn+1.)*((xn+1.)**2-xm**2)/(2.*xn+3.) ) 99 continue c now do top m=n exactly nn=2*n-1 p(n,n)=s**n do 97 i=1,nn,2 p(n,n)=p(n,n)*float(i) 97 continue y(n,n)=((-1.)**n)*sqrt( ((2*n+1)/4./pi)*fac(n-n)/fac(n+n) ) & *p(n,n)*cmplx(cos(n*phi),sin(n*phi)) c now do second to the top m=n-1 using the exact m=n, and the recursive m=n-2 found c previously xm=float(n-1) xn=float(n) y(n,n-1)=-cmplx(cos(phi),sin(phi))*y(n,n-2)* & ( xn*(xn+1.)-xm*(xm-1.) ) & /sqrt((xn+xm)*(xn-xm+1.)) y(n,n-1)=y(n,n-1)-cmplx(cos(phi),-sin(phi))*y(n,n) & *sqrt((xn-xm)*(xn+xm+1.)) y(n,n-1)=y(n,n-1)*s/2./xm/x 100 continue c now fill in -m terms do 220 n=0,nnn do 220 m=-1,-n,-1 mm=-m y(n,m)=((-1.)**mm)*cmplx(real(y(n,mm)),-aimag(y(n,mm))) 220 continue 2222 continue return end subroutine surharm(theta,phi,nnn,p,y,yphi,ytheta,yphiphi, & ythth,yphith) c computes spherical harmonics (complex) c for a value of x = cos(theta), phi=angle phi c so -1 < x < 1, P(n,m), -n < m < n, 0 < n . c uses the spherical harmonics calculated in harm real p(0:200,0:500),x,s,phi complex y(0:200,-250:250),yphith(0:200,-250:250) complex yphi(0:200,-250:250),yphiphi(0:200,-250:250) complex ytheta(0:200,-250:250),ythth(0:200,-250:250) real pqa(100),theta integer ipqa(100) pi=4.0*atan(1.0) x=cos(theta) s=sqrt(1.-x*x) c Define various partial derivatives of the y(n,m)'s using the y(n,m) c functions themselves (following Arfken) 2222 continue do 300 n=0,nnn c does n=0 separately, derivative of n=0 constant term is identically zero if(n.eq.0) then yphi(0,0)=cmplx(0.0,0.0) yphiphi(0,0)=cmplx(0.0,0.0) yphith(0,0)=cmplx(0.0,0.0) ytheta(0,0)=cmplx(0.0,0.0) ythth(0,0)=cmplx(0.0,0.0) goto 300 end if c do m=n separately and exactly, below works for all other values of m do 200 m=n-1,-n,-1 xn=float(n) xm=float(m) aaa=sqrt((2.*xn+1.)*(xn**2-xm**2)/(2.*xn-1.)) yphi(n,m)=cmplx(0.,1.)*m*y(n,m) yphiphi(n,m)=-(m**2)*y(n,m) ytheta(n,m)=xn*x*y(n,m)/s-aaa*y(n-1,m)/s ythth(n,m)=((xm*xm-xn*x*x)/s/s-xn*(xn+1.))*y(n,m)+ & aaa*x*y(n-1,m)/s/s yphith(n,m)=ytheta(n,m)*cmplx(0.,1.)*m 200 continue c now do m=n exactly m=n xn=float(n) xm=float(m) yphi(n,m)=cmplx(0.,1.)*m*y(n,m) yphiphi(n,m)=-(m**2)*y(n,m) ytheta(n,m)=xn*x*y(n,m)/s ythth(n,m)=xn*(xn*x*x-1.)*y(n,m)/s/s yphith(n,m)=ytheta(n,m)*cmplx(0.,1.)*m 300 continue return end c subroutine to make image of particle from coefficients of c spherical harmonics subroutine image(a,nnn,y,nx,ny,nz,p) real p(0:200,0:500),pi complex y(0:200,-250:250),a(0:200,-250:250),rr integer pix,edge2 pi=4.0*atan(1.0) xc=nx/2.+0.01 yc=ny/2.+0.01 zc=nz/2.+0.01 do 10 k=1,nz do 10 j=1,ny do 10 i=1,nx pix=1 23 x1=i y1=j z1=k r=sqrt( (x1-xc)**2+(y1-yc)**2+(z1-zc)**2 ) if(r.eq.0.0) then pix=3 goto 22 end if theta=acos((z1-zc)/r) phi=atan( (y1-yc)/(x1-xc) ) if((y1-yc).lt.0.0.and.(x1-xc).lt.0.0) phi=phi+pi if((y1-yc).gt.0.0.and.(x1-xc).lt.0.0) phi=phi+pi if((y1-yc).lt.0.0.and.(x1-xc).gt.0.0) phi=phi+2.*pi if(theta.ge.pi.or.theta.le.0.0) write(7+me,*) & ' error in theta, outside range & in subroutine harm',theta call harm(theta,phi,nnn,p,y) rr=cmplx(0.0,0.0) do 20 n=0,nnn do 20 m=-n,n rr=rr+a(n,m)*y(n,m) 20 continue if(r.le.real(rr)) pix=3 22 write(11+me,1) pix 1 format(i1) 10 continue return end