c computes projection of particle on z=0 plane, fits Fourier series to it, then
c computes perimeter, area, caliper dimensions in same way as slicez0-mpi-SH.f does
c for the z=0 slice
c Algorithm is to, for each value of phi  (Gaussian quadrature), compute the maximum
c in r(theta,phi)sin(theta) - this is the part that "sticks" out the most to give
c the projection. 
c **This version was made for the UConn data, and no need to save
c Fourier coefficients and perimeter points

         include 'mpif.h'

         real p(0:200,0:500),maxsize,minr,maxr,x1,y1,x2,y2,z1,z2
c aa, bb are cos and sin Fourier coefficienti
         real xg(1000),wg(1000),pi,xx1,yy1,xx,yy,aa(0:1000),bb(0:1000)
         real I11,I22,I33,I12,I13,I23,distsave(1000000,2)
         real hh,kk,rr1,rr2,meanrc,len,lineeq(1000000,6)
         real xtry,ytry,rtry,theta,phi,phipt,rprmax,xproj
         real endpts(10000,8),rf(1000),area,perim,kkurv
         real r,rphi,rphiphi,kurv(1000000),tan(1000000,3)
         real alpha,beta,gamma

c added for xcmin computation
      real xc(1000,1000,4),xcmax(1000),yk(1000),yF(100000)
      real yB(100000),phiB(100000),phiF(100000),rBB(100000),rFF(100000)
      real rpp(1000)
c note rf(j) is used for projection points around edge of 
c projection to get Fourier series

c lineeq: 1 = slope, 2 = intercept. 3,4 - x,y of one end of segment, 5,6 = other end
         real slope,bint
         complex y(0:200,-250:250),a(0:200,-250:250)
         complex b(0:200,-250:250),a3(0:200,-250:250)
         complex yphi(0:200,-250:250),sa
         complex ytheta(0:200,-250:250)
         complex ythth(0:200,-250:250),yphith(0:200,-250:250)
         complex f1,yphiphi(0:200,-250:250),r1
         complex ze,zf,zg,zl,zm,zn,zq,zh,zh2,zk,zq1
         complex ddd
         complex thx,thy,thz,phix,phiy,phiz
         complex nxphi,nxth,nyphi,nyth,nzphi,nzth
         complex testhh,testkk,inttest,rr
         integer pix(200,200)
         integer ierr,npes,me,Ntot,first,last,ios,first1,first2
         integer*4 iseedc(2)

c USER
         character*23 namef
         character*26 filename
         character*30 filename2
         character*50 namef13
         character*31 namefproj
         character*27 namefproj2

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=================================================================

         filename=''

c USER
         open(unit=7,file='ASTM-AMPM-2-proj-SH-data.dat')

c USER
c set the number of anm files
         Ntot=14581

         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 routine mpi time to know how long each processor takes to run
         tc0=MPI_Wtime()
         tc1=MPI_Wtime()
         tar=tc1-tc0
         tc0=MPI_Wtime()

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)
         end if

c root processor sends a 2 vector data to all other processors
c the vector data refers to xg and wg read from Gaussian 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=======================================================================

c program reads from this geom file, gets anm filename and other parameters
         open(unit=29,file='ASTM-AMPM2-geom-len.dat')

c USER
c  creation of output files that depend on the process rank
c this will eventually be the file that has basic information
c on each projection, but not the list of 2D Fourier coefficients
         write(namef(16:19),'(I4)') 1000+me
         namef(1:16)='ASTM-AMPM2-proj-'
         namef(20:23)='.dat'

         open(unit=12,file=namef)

c calculation of the number of files that each processor must read           
        last= Ntot/npes
        first1= Ntot/npes 
        first2=me*first1
        rem = Ntot-npes*first1
        if (rem.ne.0) then
          if(me.ge.(npes-rem)) then
          last=last+1
          first2=first2+(me-(npes-rem))
          end if
        end if

c set the first file that each processor must read
c have each processor read down to where it should start
c except me = 0

         if(me.gt.0) then

         do 155 ijk=1,first2
         read(29,1888) filename
c USER
1888     format(a26)
155      continue
         end if

        iseedc(1)=(me+13)*193
        call random_seed(put=iseedc)

c loop in particles
         do 3800 lmn=1,last
c read information about particles from geom file
c USER
         read(29,299) filename,nnn
299      format(a26,148x,i2)
c         namefproj(1:14)=filename
c         namefproj(11:31)='-Fourier-coeff-z0.txt'
c         open(unit=21,file=namefproj)
c         namefproj2(1:14)=filename
c         namefproj2(11:27)='-2D-points-z0.txt'
c         open(unit=22,file=namefproj2)
c USER
         filename2(5:30)=filename
         filename2(1:4)='anm/'

         open(unit=19,file=filename2)
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 only read in enough for this particle
         do 390 n=0,nnn
         do 390 m=n,-n,-1
         read(19,*) ii,jj,a1,a2
         a(n,m)=cmplx(a1,a2)
c aa is original, save through do 3888 loop
         a3(n,m)=cmplx(a1,a2)
390      continue
         close(19)

c loop in x,y,z projections, so three for each particle
         do 3888 iijjkk=1,3
c choose x, or y, or z axis for projection, randomly
c z all angles zero
c x beta = pi/2
c y alpha=pi/2  beta = pi/2

c        call random_number(xrot1)
c        call random_number(xrot2)

c  Rotate coefficients, a(n,m) into b(n,m)
c always do projection in z direction, but first rotate 
c particle so we get the x and y directions

c choose all zero angles to do projection in z direction
      if(iijjkk.eq.3) then
         alpha=0.0
         beta=0.0
         gamma=0.0
       end if

c rotate y into z to get projection in y direction
c y is now z-axis, so get projection in y direction
      if(iijjkk.eq.2) then
         alpha=-0.5*pi
         beta=-0.5*pi
         gamma=0.0
       end if

c rotate x into z to get projection in x direction
c x is now in z-direction, so get projection in z direction
      if(iijjkk.eq.1) then
         alpha=pi
         beta=-0.5*pi
         gamma=0.0
       end if

         do 520 n=0,nnn
         do 520 m=-n,n
         b(n,m)=cmplx(0.,0.)
          do 510 mp=-n,n
          call drot(pi,alpha,beta,gamma,n,m,mp,ddd)
          b(n,m)=b(n,m)+a3(n,mp)*ddd
510       continue
520       continue

c now copy b(n,m) back into a(n,m), since all code below uses a(n,m)

         do 521 n=0,nnn
         do 521 m=-n,n
          a(n,m)=b(n,m)
521       continue

c compute projection points in z=0 plane
         ntheta=200
         do 830 j=1,ng
         rprmax=-10000.0
         phi=pi*xg(j)+pi
         do 820 i=1,ntheta
         theta=pi*float(i)/float(ntheta+1)
         call harm(theta,phi,nnn,p,y)

         rr=a(0,0)*y(0,0)
c  derivatives of (0,0) term are zero, since it is a constant
         do 22 n=1,nnn
         do 22 m=n,-n,-1
         rr=rr+a(n,m)*y(n,m)
22       continue
         
         xproj=real(rr)*sin(theta)
         if(rprmax.lt.xproj) rprmax=xproj

820      continue
c point for projection in z=0 plane for phi(j)
c rf(j) is the 2D set of points on the projection surface
         rf(j)=rprmax
c         write(22,833) rprmax*cos(phi),rprmax*sin(phi)
c833      format(2f20.10)
830      continue
         close(22)

c now fit Fourier series to these values
         
c  nnnF is how many sin and cos terms to use in series (n limit)
c bb(0)=0.0
         nnnF=nnn
         do 391 n=0,nnnF
         aa(n)=0.0
         bb(n)=0.0
391      continue

         do 831 j=1,ng
         phi=pi*xg(j)+pi
         aa(0)=aa(0)+rf(j)*wg(j)
         do 400 n=1,nnnF
         aa(n)=aa(n)+rf(j)*wg(j)*cos(n*phi)
         bb(n)=bb(n)+rf(j)*wg(j)*sin(n*phi)
400      continue

831      continue

         do 832 j=1,ng
         phi=pi*xg(j)+pi
         r=aa(0)*0.5
         do 402 n=1,nnnF
         r=r+aa(n)*cos(n*phi)+bb(n)*sin(n*phi)
402      continue
         
832      continue

c         do 404 n=0,nnnF
c         write(21,219) n,aa(n),bb(n)
c219      format(i5,2f20.12)
c404      continue
c         close(21)

c so at this point, aa(n),bb(n), n=0,nnnF are the Fourier coefficients
         area=0.0
         perim=0.0

         factor=pi
c redo x and y extents, since projection is different from a slice
c projection should always  have the same extents in x and y, if
c projected in z=0 plane
         x1=10000.
         y1=10000.
         x2=0.0
         y2=0.0
         do 430 j=1,ng
         phi=pi*xg(j)+pi

         r=aa(0)*0.5
         rphi=0.0
         rphiphi=0.0
c  derivatives of(0,0) term is zero, since it is a constant
         do 62 n=1,nnnF
         r=r+aa(n)*cos(n*phi)+bb(n)*sin(n*phi)
         rphi=rphi+n*(-aa(n)*sin(n*phi)+bb(n)*cos(n*phi))
         rphiphi=rphiphi-n*n*(aa(n)*cos(n*phi)+bb(n)*sin(n*phi))
62       continue
         xx1=r*cos(phi)
         yy1=r*sin(phi)
         if(x1.gt.xx1) x1=xx1*1.00001
         if(y1.gt.yy1) y1=yy1*1.00001
         if(x2.lt.xx1) x2=xx1*1.00001
         if(y2.lt.yy1) y2=yy1*1.00001
         
         perim=perim+wg(j)*sqrt(rphi**2+r**2)
         area=area+0.5*wg(j)*r*r
430      continue
         dmin=0.25*min(x2-x1,y2-y1)
         area=area*factor
         perim=perim*factor
         if(me.eq.0) then
         write(7,*) filename2,area,perim
         write(7,*) sqrt(4.*area/pi),perim/pi
         end if


c go around object, use Nj angles, remember those with positive curvature
         Nj=5000
         jcurv=0
c initialize lineeq line segment end points so I know which have been used
         do 489 i=1,1000000
         lineeq(i,3)=-10000.
         lineeq(i,4)=-10000.
         lineeq(i,5)=-10000.
         lineeq(i,6)=-10000.
489      continue

         do 530 j=1,Nj
         phi=(j-1)*2.*pi/float(Nj)

         r=aa(0)*0.5
         rphi=0.0
         rphiphi=0.
c  derivatives of (0,0) term is zero, since it is a constant
         do 32 n=1,nnnF
         r=r+aa(n)*cos(n*phi)+bb(n)*sin(n*phi)
         rphi=rphi+n*(-aa(n)*sin(n*phi)+bb(n)*cos(n*phi))
         rphiphi=rphiphi-n*n*(aa(n)*cos(n*phi)+bb(n)*sin(n*phi))
32       continue
         xx1=r*cos(phi)
         yy1=r*sin(phi)
         if(xx1.lt.x1) x1=1.00001*xx1
         if(xx1.gt.x2) x2=1.00001*xx1
         if(yy1.lt.y1) y1=1.00001*yy1
         if(yy1.gt.y2) y2=1.00001*yy1
         
         kkurv=(r*r+2.*(rphi)**2-r*rphiphi)
         kkurv=kkurv/(r*r+rphi**2)**(1.5)

         if(kkurv.ge.0.0) then
          jcurv=jcurv+1
          kurv(jcurv)=kkurv
c work out tangent unit vector
          tan(jcurv,1)=(rphi*cos(phi)-r*sin(phi))/sqrt(rphi**2+r**2)
          tan(jcurv,2)=(rphi*sin(phi)+r*cos(phi))/sqrt(rphi**2+r**2)
          tan(jcurv,3)=tan(jcurv,2)/tan(jcurv,1)
c work out slope and y-intercept of tangent line
          slope=tan(jcurv,3)
          bint=r*sin(phi)-tan(jcurv,3)*r*cos(phi)
          lineeq(jcurv,1)=slope
          lineeq(jcurv,2)=bint

c check intersection at top or bottom, then left or right
c this assumes that slope is never vertical or zero

c this is the first if statement, no end points have yet been found, so if one 
c is found here, it goes in (3,4)
          xx=(y1-bint)/slope
          if(xx.ge.x1.and.xx.le.x2) then
           lineeq(jcurv,3)=xx
           lineeq(jcurv,4)=y1
          end if

          xx=(y2-bint)/slope
          if(xx.ge.x1.and.xx.le.x2) then
           if(lineeq(jcurv,3).eq.-10000.) then
            lineeq(jcurv,3)=xx
            lineeq(jcurv,4)=y2
            goto 1200
           end if
            lineeq(jcurv,5)=xx
            lineeq(jcurv,6)=y2
          end if

1200      yy=slope*x1+bint
          if(yy.ge.y1.and.yy.le.y2) then
           if(lineeq(jcurv,3).eq.-10000.) then
            lineeq(jcurv,3)=x1
            lineeq(jcurv,4)=yy
            goto 1300
           end if
           if(lineeq(jcurv,5).eq.-10000.) then
            lineeq(jcurv,5)=x1
            lineeq(jcurv,6)=yy
           end if
          end if

1300      yy=slope*x2+bint
          if(yy.ge.y1.and.yy.le.y2) then
           if(lineeq(jcurv,3).eq.-10000.) then
            lineeq(jcurv,3)=x2
            lineeq(jcurv,4)=yy
           end if
           if(lineeq(jcurv,5).eq.-10000.) then
            lineeq(jcurv,5)=x2
            lineeq(jcurv,6)=yy
           end if
          end if

c now endpoints of line segment are established, determine if points on line
c segment are inside particle or not. If not, then keep this point. If yes,
c then throw away this point by letting jcurv=jcurv-1
c does not check actual endpoints, only points in-between
c assumes origin is at center of particle (which is how the extents, anm, etc. are defined

         nover=100
         do 475 it=1,nover-1
         t=float(it)/float(nover)
         xtry=lineeq(jcurv,3)*(1.-t)+t*lineeq(jcurv,5)
         ytry=lineeq(jcurv,4)*(1.-t)+t*lineeq(jcurv,6)
         rtry=sqrt(xtry**2+ytry**2)
       phipt=atan( ytry/xtry )
       if(ytry.lt.0.0.and.xtry.lt.0.0) phipt=phipt+pi
       if(ytry.gt.0.0.and.xtry.lt.0.0) phipt=phipt+pi
       if(ytry.lt.0.0.and.xtry.gt.0.0) phipt=phipt+2.*pi
       r=0.5*aa(0)
       do 20 n=1,nnnF
       r=r+aa(n)*cos(n*phipt)+bb(n)*sin(n*phipt)
20     continue
        if(rtry.le.r) then
         lineeq(jcurv,3)=-10000.
         lineeq(jcurv,4)=-10000.
         lineeq(jcurv,5)=-10000.
         lineeq(jcurv,6)=-10000.
         jcurv=jcurv-1
        goto 530
        end if
475      continue

c this is the end of curvature > 0 if then statement
         end if

530      continue

c if program has made it to here, then jcurv is the total number of slopes at positive
c curvature points, and line segment to bounding box does not overlap the particle
c now compare pairs of these to find calipers, take maximum of pairs for each slope considered
c to get the caliper distance for that slope of tangent lines
         if(me.eq.0) write(7,*) ' number of good points = ',jcurv
         npp=0
         distmax=-10.0
         distmin=10000.0
         do 610 j=1,jcurv-1
         np=0
c comparing slopes, if difference is more than 1 %, go to next pair
         do 600 jj=j+1,jcurv
         diff=abs( (lineeq(j,1)-lineeq(jj,1))/lineeq(j,1))
         if(diff.gt.0.01) goto 600
         dist=abs(lineeq(j,2)-lineeq(jj,2))/sqrt(1.+lineeq(j,1)**2)
         if(dist.lt.dmin) goto 600
         if(distmax.lt.dist) then
         np=np+1
         distmax=dist
         j1=j
         j2=jj
         end if
         if(distmin.gt.dist) then
         np=np+1
         distmin=dist
         jj1=j
         jj2=jj
         end if
600      continue

610      continue
         Fermax=distmax
         Fermin=distmin

c now compute xcmin
c xc(m,k,2) m=phi0s, k=y's, 1=phiF, 2=phiB (x - coordinates)
c xcmax(m) = max over k of abs(xc(m,k,1)-xc(m,k,2))
c xcmin = min over m of xcmax(m)

c loop in rotation angle phi0 (index m)
c only need to do 0-pi, since particle up is the same as particle down
c this is different orientations of the projection
c Note - think of this as: rprime(phi)=r(phi-phi0)
c so when I compute the corresponding y value, I use phi, not phi-phi0
       Mphi0=200
       do 4180 m=1,Mphi0
c rotation angle
       phi0=(m-1)*pi/float(Mphi0)

c compute ymin and ymax and angles at these values for given phi0
       ymin=10000.
       ymax=0.0
c Nex = number of angles used to calculate ymin and ymax
c phi is the true angle, phi-phi0 gives the rotated version
       Nex=1000
       do 4130 j=1,Nex
       phi=(j-1)*2.*pi/float(Nex)
c       if(phi.lt.0.0) phi=phi+2.*pi
c       if(phi.gt.2.*pi) phi=phi-2.*pi
       r=aa(0)*0.5
       ppp=phi-phi0
       do 612 n=1,nnnF
       r=r+aa(n)*cos(n*ppp)+bb(n)*sin(n*ppp)
612    continue
       yy1=r*sin(phi)
       if(ymin.gt.yy1) then
       ymin=yy1
       phimin=phi
       end if
       if(ymax.lt.yy1) then
       ymax=yy1
       phimax=phi
       end if
4130   continue
c       if(phimax.lt.0.0) phimax=phimax+2.*pi
c       if(phimax.gt.2.*pi) phimax=phimax-2.*pi
c       if(phimin.lt.0.0) phimin=phimin+2.*pi
c       if(phimin.gt.2.*pi) phimin=phimin-2.*pi
c       if(me.eq.0) write(33,*) filename,phi0,phimin,ymin,phimax,ymax

c  ymin,phimin and ymax,phimax are now computed for m (phi0) orientation
c lay out yk(k) array for KK-1 points
       KK=200
       do 4131 k=1,KK-1
       yk(k)=ymin+k*(ymax-ymin)/float(KK)
       if(yk(k).eq.0.0) yk(k)=yk(k)+0.1*(ymax-ymin)/float(KK)
c       if(me.eq.0) write(33,*) filename,area,k,yk(k)
4131   continue

c now determine angles needed to match set of yk(k)'s
c start with angle phimin, and go both ways until reaching phimax
      NNphi=10000
c F = forward (counter clockwise) B = backwards (clockwise)
c increment in angle for backward direction
      phiincB=abs(phimax-phimin)/float(NNphi)
c increment in angle for forward direction
      phiincF=(2.*pi-abs(phimax-phimin))/float(NNphi)
      kF=0
      kB=0
      do 4150 nphi=1,NNphi-1
      phiF(nphi)=phimin+nphi*phiincF
      phiB(nphi)=phimin-nphi*phiincB
c      if(phiB(nphi).lt.0.0) phiB(nphi)=phiB(nphi)+2.*pi
c      if(phiB(nphi).gt.2.*pi) phiB(nphi)=phiB(nphi)-2.*pi
c      if(phiF(nphi).lt.0.0) phiF(nphi)=phiF(nphi)+2.*pi
c      if(phiF(nphi).gt.2.*pi) phiF(nphi)=phiF(nphi)-2.*pi
      rFF(nphi)=aa(0)*0.5
      rBB(nphi)=aa(0)*0.5
c need to rotate angles when calculating rFF and rBB
      aaa=phiF(nphi)-phi0
      bbb=phiB(nphi)-phi0
c      if(aaa.lt.0.0) aaa=aaa+2.*pi
c      if(aaa.gt.2.*pi) aaa=aaa-2.*pi
c      if(bbb.lt.0.0) bbb=bbb+2.*pi
c      if(bbb.gt.2.*pi) bbb=bbb-2.*pi
      do 613 n=1,nnnF
      rFF(nphi)=rFF(nphi)+aa(n)*cos(n*aaa)+bb(n)*sin(n*aaa)
      rBB(nphi)=rBB(nphi)+aa(n)*cos(n*bbb)+bb(n)*sin(n*bbb)
613   continue
c use unrotated angle when computing y coordinate
      yF(nphi)=rFF(nphi)*sin(phiF(nphi))
      yB(nphi)=rBB(nphi)*sin(phiB(nphi))

      if(yF(nphi).lt.yk(kF+1)) then
       if(yB(nphi).lt.yk(kB+1)) then
       goto 4150
       end if
      end if

      if(kF.eq.(KK-1)) goto 4114
      if(yF(nphi).gt.yk(kF+1)) then
      xm=(yF(nphi)-yF(nphi-1))/(phiF(nphi)-phiF(nphi-1))
      phiFsave=phiF(nphi-1)+(yk(kF+1)-yF(nphi-1))/xm
c compute r(phiFsave)
      aaa=phiFsave-phi0
c      if(aaa.lt.0.0) aaa=aaa+2.*pi
c      if(aaa.gt.2.*pi) aaa=aaa-2.*pi
      rrF=aa(0)*0.5
      do 614 n=1,nnnF
      rrF=rrF+aa(n)*cos(n*aaa)+bb(n)*sin(n*aaa)
614   continue
      kF=kF+1
      xc(m,kF,1)=rrF*cos(phiFsave)
      xc(m,kF,3)=rrF*sin(phiFsave)
c      if(me.eq.0) then
c        write(33,*) 'F ',filename,area,kF,xc(m,kF,1),xc(m,kF,3)
c      end if
      end if

4114  if(kB.eq.(KK-1)) goto 4115
      if(yB(nphi).gt.yk(kB+1)) then
      xm=(yB(nphi)-yB(nphi-1))/(phiB(nphi)-phiB(nphi-1))
      phiBsave=phiB(nphi-1)+(yk(kB+1)-yB(nphi-1))/xm
c compute r(phiBsave)
      bbb=phiBsave-phi0
c      if(bbb.lt.0.0) bbb=bbb+2.*pi
c      if(bbb.gt.2.*pi) bbb=bbb-2.*pi
      rrB=aa(0)*0.5
      do 617 n=1,nnnF
      rrB=rrB+aa(n)*cos(n*bbb)+bb(n)*sin(n*bbb)
617   continue
      kB=kB+1
      xc(m,kB,2)=rrB*cos(phiBsave)
      xc(m,kB,4)=rrB*sin(phiBsave)
c      if(me.eq.0)then
c        write(33,*) 'B ',filename,area,kB,xc(m,kB,2),xc(m,kB,4)
c      end if
      end if
4115  if(kB.eq.(KK-1).and.kF.eq.(KK-1)) goto 5555

4150  continue
c5555  if(me.eq.0) write(33,*) 'kB kF after 4150 loop ',kB,kF
5555  continue

c now find maximum over kB and kF to determine xcmax(m)
      xcmax(m)=0.0
      do 4155 k=1,KK-1
      dd=abs(xc(m,k,1)-xc(m,k,2))
      if(xcmax(m).lt.dd) xcmax(m)=dd
c      write(32,*) filename,area,xc(m,k,1),xc(m,k,3),xc(m,k,2),
c     &      xc(m,k,4),xcmax(m) 
4155  continue
c go to another rotation angle
4180  continue

      xcmin=10000.0
c xcm = max of xcmax
      xcm=0.0
      do 4190 m=1,Mphi0
      if(xcmin.gt.xcmax(m)) xcmin=xcmax(m)
      if(xcm.lt.xcmax(m)) xcm=xcmax(m)
4190  continue

c now xcmin has been calculated, write out all results

      write(12,122) filename,area,perim,Fermin,Fermax,xcmin,xcm,nnn

122   format(a26,6f18.6,i4)
c end of x,y,z projection loopfor a given particle
3888     continue

3800     continue
c end of particle loop
      close(12)
      close(29)

c calculation of time that each processor has used
      tc1=MPI_Wtime()
      sec = (tc1-tc0) - tar
      print *,'proc = ',me,' finished at ',sec,' seconds' 

c mpi barrier used here to make sure that all processors should finish together
      call mpi_barrier(MPI_COMM_WORLD,ierr)

c===================================================================
c now concatenate the results files
      if(me.eq.0) then

c USER
      do 2201 i=0,npes-1

           write(namef13(20:23),'(I4)') 1000+i
           namef13(1:20)='cat ASTM-AMPM2-proj-'
           namef13(24:50)='.dat >> ASTM-AMPM2-proj.dat'

c runs cat statement
         call system(namef13)

2201       continue

           end if

         call mpi_finalize(ierr)

         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)
c       ytheta(n,m)=xn*x*y(n,m)/s-aaa*y(n-1,m)/s
c       ythth(n,m)=((xm*xm-xn*x*x)/s/s-xn*(xn+1.))*y(n,m)+
c     &             aaa*x*y(n-1,m)/s/s
c       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)
c       ytheta(n,m)=xn*x*y(n,m)/s
c       ythth(n,m)=xn*(xn*x*x-1.)*y(n,m)/s/s
c       yphith(n,m)=ytheta(n,m)*cmplx(0.,1.)*m
300    continue
     
       return
       end
      subroutine drot(pi,alpha,beta,gamma,n,m,mp,ddd)
      real alpha,beta,gamma,pi
      complex ddd

      cosbeta=cos(beta/2.)
      sinbeta=sin(beta/2.)
c  allow for correct limiting behavior of (0)^0 = 1
      if(cosbeta.eq.0.0) then
       beta=beta+1.0e-10
       cosbeta=cos(beta/2.)
      end if
      if(sinbeta.eq.0.0) then
       beta=beta+1.0e-10
       sinbeta=sin(beta/2.)
      end if
    
      ddd=sqrt(fac(n+mp)*fac(n-mp)/fac(n+m)/fac(n-m))
      klow=max(0,m-mp)
      khigh=min(n-mp,n+m)
      
      total=0.0
      do 100 k=klow,khigh
      abc=((-1)**(k+mp-m))*( fac(n+m)/fac(k)/fac(n+m-k) )*
     &     ( fac(n-m)/fac(n-mp-k)/fac(-m+mp+k) )
      total=total+abc*(cosbeta**(2*n+m-mp-2*k))*
     &      (sinbeta**(2*k+mp-m))
100   continue
      ddd=ddd*total*cmplx(cos(mp*alpha),-sin(mp*alpha))*
     &     cmplx(cos(m*gamma),-sin(m*gamma))
      
      return
      end