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