c Computes L, W, and T for nonSH files, from the voxel representation c new aspect - sx is now in the *-part-*.dat file along with all c the voxels real pi,sx real xn,yn,zn,alphaL,betaL real wwwx,wwwy,wwwz,alphaW,betaW real tttx,ttty,tttz,alphaT,betaT integer save(10000000,4) c needs to be changed with filename length changes character*33 filename character*39 filename2 integer ierr,npes,me,Ntot,first,last,ios,first1 integer nnn1 integer n2,ind,ss,ss2 integer rem,nnsave,iii real length,width,thick,leth,with,thth c USER c set the number of reject files Ntot=2390 c USER c voxel length scale (micrometers) is sx pi=4.0*atan(1.0) nnsave=10000000 iii=4 open(unit=29,file='nonSH.lis') c USER c replace characters before "-reject" c open(unit=12,file='Num-1-Praxair-nonSH-len.dat') open(unit=24,file='ASTM-AMPM2-L-nonSH-unit-vector.dat') open(unit=26,file='ASTM-AMPM2-W-nonSH-unit-vector.dat') open(unit=28,file='ASTM-AMPM2-T-nonSH-unit-vector.dat') do 3800 lmn=1,Ntot read(29,2999) filename c USER 2999 format(a33) filename2(1:6)='nonSH/' filename2(7:39)=filename open(unit=19,file=filename2) c Read in voxel coordinates (middle of voxel) c first read in number of voxels and scale (um) read(19,1003) isave,sx 1003 format(i10,f9.6) vol=isave*sx*sx*sx do 390 n=1,isave read(19,*) save(n,1),save(n,2),save(n,3),save(n,4) 390 continue close(19) call lwt(sx,isave,save,length,width,thick,dotlw,dotlt,dotwt, & nnsave,iii, & xn,yn,zn,alphaL,betaL, & wwwx,wwwy,wwwz,alphaW,betaW, & tttx,ttty,tttz,alphaT,betaT) write(12,399) filename,vol,length,width,thick,dotlw,dotlt,dotwt, & xn,yn,zn,alphaL,betaL, & wwwx,wwwy,wwwz,alphaW,betaW, & tttx,ttty,tttz,alphaT,betaT 399 format(a33,f18.8,3f14.8,3f7.3,15f13.8) call flush(12) c lwt subroutine done 3800 continue 1001 close(12) close(29) end c ====================================================================================================== subroutine lwt(sx,isave,save,length,width,thick,dotlw,dotlt,dotwt, & nnsave,iii, & xn,yn,zn,alphaL,betaL, & wwwx,wwwy,wwwz,alphaW,betaW, & tttx,ttty,tttz,alphaT,betaT) c Computes length (maximum distance between two surface points on particle), then maximum c distance (width) between two points that is perpendicular to length, and then thickness, c which is the maximum distance perpendicular to both length and thickness. c Uses voxel representation of particle real length,width,thick,pi,sx,dotlw,dotlt,dotwt real xn,yn,zn,alphaL,betaL real wwwx,wwwy,wwwz,alphaW,betaW real tttx,ttty,tttz,alphaT,betaT integer nnsave,iii integer save(nnsave,iii) 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) pi=4.0*atan(1.0) rmax=0.0 do 5555 i2=1,isave-1 do 5555 j2=i2+1,isave if(save(i2,4).eq.0.or.save(j2,4).eq.0) goto 5555 c Compute length xx1=save(i2,1) yy1=save(i2,2) zz1=save(i2,3) xx2=save(j2,1) yy2=save(j2,2) zz2=save(j2,3) rm=sqrt((xx1-xx2)**2+(yy1-yy2)**2+(zz1-zz2)**2) if(rmax.lt.rm) then rmax=rm ns1=i2 ns2=j2 end if 5555 continue length=rmax c (xn,yn,zn) is unit vector along length c from 1 to 2 xn=(save(ns2,1)-save(ns1,1))/length yn=(save(ns2,2)-save(ns1,2))/length zn=(save(ns2,3)-save(ns1,3))/length rmax=0.0 c now compute width do 5557 i2=1,isave-1 do 5557 j2=i2+1,isave if(save(i2,4).eq.0.or.save(j2,4).eq.0) goto 5557 c Compute length xx1=save(i2,1) yy1=save(i2,2) zz1=save(i2,3) xx2=save(j2,1) yy2=save(j2,2) zz2=save(j2,3) 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 c only widths within acos(0.05) or > 87 degrees are considered if(abs(dot).le.0.05) then if(rmax.lt.rm) then rmax=rm ns3=i2 ns4=j2 wwx=wx wwy=wy wwz=wz end if end if 5557 continue width=rmax ddot=wwx*xn+wwy*yn+wwz*zn dot=acos(ddot) c set final width unit vector and angles to use for thickness wwwx=wwx wwwy=wwy wwwz=wwz rmax=0.0 c thickness: (xn,yn,zn) and (wwwx,wwwy,wwwz) are length and width unit vectors do 5559 i2=1,isave-1 do 5559 j2=i2+1,isave if(save(i2,4).eq.0.or.save(j2,4).eq.0) goto 5559 c Compute thickness xx1=save(i2,1) yy1=save(i2,2) zz1=save(i2,3) xx2=save(j2,1) yy2=save(j2,2) zz2=save(j2,3) 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 if(abs(dotl).le.0.05.and.abs(dotw).le.0.05) then c if(abs(dotl).le.0.25.and.abs(dotw).le.0.25) then if(abs(dotl).le.0.15.and.abs(dotw).le.0.15) then if(rmax.lt.rm) then rmax=rm ns5=i2 ns6=j2 ttx=tx tty=ty ttz=tz end if end if 5559 continue thick=rmax c set final thickness unit vector tttx=ttx ttty=tty tttz=ttz thick=thick*sx width=width*sx length=length*sx c dot products for LW, LT, WT dotlw=(180./pi)*acos(wwwx*xn+wwwy*yn+wwwz*zn)-90. dotlt=(180./pi)*acos(tttx*xn+ttty*yn+tttz*zn)-90. dotwt=(180./pi)*acos(wwwx*tttx+wwwy*ttty+wwwz*tttz)-90. 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