22 subroutine clustering(xl,yl,zl,n,xclust,yclust,zclust,fclust,rms, &
61 integer :: n,i,j,l,nclust(maxpart),numb(ncluster),ncl
62 real :: xl(n),yl(n),zl(n),xclust(ncluster),yclust(ncluster),x,y,z
63 real :: zclust(ncluster),
distance2,distances,distancemin,rms,rmsold
64 real :: xav(ncluster),yav(ncluster),zav(ncluster),fclust(ncluster)
65 real :: rmsclust(ncluster)
70 if (n.lt.ncluster)
return
88 xclust(j)=xl(j*n/ncluster)
89 yclust(j)=yl(j*n/ncluster)
106 distances=
distance2(yl(i),xl(i),yclust(j),xclust(j))
107 if (distances.lt.distancemin)
then
108 distancemin=distances
130 numb(nclust(i))=numb(nclust(i))+1
132 yclust(nclust(i)),xclust(nclust(i)))
138 rms=rms+distances*distances
139 rmsclust(nclust(i))=rmsclust(nclust(i))+distances*distances
144 x = cos(yl(i))*sin(xl(i))
145 y = -1.*cos(yl(i))*cos(xl(i))
147 xav(nclust(i))=xav(nclust(i))+x
148 yav(nclust(i))=yav(nclust(i))+y
149 zav(nclust(i))=zav(nclust(i))+z
152 rms=sqrt(rms/
real(n))
159 if (numb(j).gt.0)
then
160 rmsclust(j)=sqrt(rmsclust(j)/
real(numb(j)))
161 xav(j)=xav(j)/
real(numb(j))
162 yav(j)=yav(j)/
real(numb(j))
163 zav(j)=zav(j)/
real(numb(j))
168 xclust(j)=atan2(xav(j),-1.*yav(j))
169 yclust(j)=atan2(zav(j),sqrt(xav(j)*xav(j)+yav(j)*yav(j)))
177 if ((l.gt.1).and.(abs(rms-rmsold)/rmsold.lt.0.005)) goto 99
190 zclust(nclust(i))=zclust(nclust(i))+zl(i)
194 xclust(j)=xclust(j)/pi180
195 yclust(j)=yclust(j)/pi180
196 if (numb(j).gt.0) zclust(j)=zclust(j)/
real(numb(j))
197 fclust(j)=100.*
real(numb(j))/
real(n)
205 zdist=zl(i)-zclust(nclust(i))
206 zrms=zrms+zdist*zdist
208 if (zrms.gt.0.) zrms=sqrt(zrms/
real(n))
real function distance2(rlat1, rlon1, rlat2, rlon2)
subroutine clustering(xl, yl, zl, n, xclust, yclust, zclust, fclust, rms, rmsclust, zrms)