source: PSPA/parmelaPSPA/trunk/emit90dp.f @ 411

       subroutine EMIT90DP (nt90,n,dpx,dpy,dpa,dpb,e90,e100,e)
c-----------------------------------------------------------------------
c      **name**        EMIT90          filename is EMIT90P.for
c
c      **description** This is a patch file for routine emit90 in PARMELA.
c      EMIT90 calculates the 90% and 100% particle-enclosing ellipse areas.
c      This version implements a speedier algorithm for finding the 90th
c      percentile particle.  The z-array is destroyed.
c
c      Old routine: N*N/10 points scanned.
c      Present version: 4-6*N points scanned.
c
c      Further speedup is possible with an extra storage array.
c
c      **written**     1.00    07/01/86                        rajeev rohatgi
c
c      **revised**     1.01    07/01/86        working         rajeev rohatgi
c
c----------------------------------------------------------------------------
       implicit double precision (e)
c
       include 'param_sz.h'
       include 'ucom.h'
c
       double precision dpa,dpb,a,b
       double precision dpx(n),dpy(n)
       double precision dpxb,dpyb,xb,yb
       double precision x(imaa),y(imaa),z(imaa)
       double precision xi,yi,zi,s0,s1,s2,g,atmp,sk
c--------------------------------------------------------------------------
c*
       if (n.le.1) return
       do 11 i=1,n
       x(i)=dpx(i)
       y(i)=dpy(i)
 11    continue
C     FIND ELLIPSE SHAPE PARAMETERS
       call elipsdp (nt90,n,dpx,dpy,dpa,dpb,dpxb,dpyb,e)
       if(e.eq.0.0) return
c      if you want 80 percent change 0.9 in 0.8 and in format 330 of outlal
       n90     = n*0.9+0.5
       a=dpa
       b=dpb
       xb=dpxb
       yb=dpyb
       if(b.eq.0.0) return
       g       = (1.+a*a)/b
       atmp    = 2.*a
       s0      = 0.
       s1      = 0.
C     CALCULATE AN AREA ASSOCIATED WITH EACH POINT
       do 10 i=1,n
       xi      = x(i)-xb
       yi      = y(i)-yb
       zi      = xi*(g*xi+atmp*yi)+b*yi*yi
       x(i)    = xi
       y(i)    = yi
       z(i)    = zi
       s0      = s0+zi
       s1      = dmax1(zi,s1)
 10    continue
       s0      = s0/n
       nk      = n
       m       = n+1-n90
C     LOOP TO SCAN AND PARTITION PARTICLES WRT MEAN
 20    continue
       kk      = 0
C     COUNT PARTICLES ABOVE MEAN
       do 29 i=1,nk
       if (z(i).gt.s0) kk = kk+1
 29    continue
c-------------------------------------------------------------------------
       if(kk.eq.0) then
        if(nt90.eq.1) then
       write(nemit,*) 
     * ' kk = 0 for (x,xp) in emit 90 no calculation of e90 and e100'
        else
       write(nemit,*) 
     * ' kk = 0 for (y,yp) in emit 90 no calculation of e90 and e100'
        endif
       e90 = 0.
       e100 = 0.
       else
       if (kk.ge.m) then
         if (kk.eq.m) then
C                                      TAKE BOTTOM-MOST OF THE 'ABOVE' SET
           s2  = s1
           do 39 i=1,nk
           if (z(i).gt.s0) s2 = dmin1(s2,z(i))
 39        continue
         else
C                                      TAKE THE 'ABOVE' SET
           sk  = 0.
           j   = 1
           do 49 i=1,nk
           zi  = z(i)
           if (zi.le.s0) go to 49
           z(j) = zi
           sk  = sk+zi
           j   = j+1
 49        continue
           nk  = kk
           s0  = sk/nk
C                                      LOOP BACK TO PARTITION REDUCED SET
           go to 20
         end if
       else
         m     = m-kk
         if (m.eq.1) then
C                                      TAKE TOP-MOST OF 'BELOW' SET
           s2  = 0.
           do 59 i=1,nk
           if (z(i).le.s0) s2 = dmax1(s2,z(i))
 59        continue
         else
C                                      TAKE 'BELOW' SET
           sk  = 0.
           j   = 1
           do 69 i=1,nk
           zi  = z(i)
           if (zi.gt.s0) go to 69
           z(j) = zi
           sk  = sk+zi
           j   = j+1
 69        continue
           nk  = nk-kk
           s0  = sk/nk
C                                      LOOP BACK TO PARTITION REDUCED SET
           go to 20
         end if
       end if
       e90     = s2
       e100    = s1
       endif
c-----------------------------------------------------------------------
       return
       end