C C This code is distributed under the terms and conditions of the C CCP4 licence agreement as `Part ii)' software. See the conditions C in the CCP4 manual for a copyright statement. C C C*********************** C Program FSEARCH ** C ********************** c Revision C Original copy:f3.f from crystal C Purpose: 6D search C ******INITIALISE******* PROGRAM FSEARCH PARAMETER(MAXREF=120000) PARAMETER(MAXPRO=200000) PARAMETER(MAX3D=5324800) parameter(npeak=200) COMPLEX DATA(MAX3D),FC(MAXREF,64) real ed(MAX3D) COMPLEX FCOMP(MAXREF),shift REAL RES(MAXREF),FCSUM(MAXREF),FF(MAXREF) INTEGER NN(3),CHKOV INTEGER*2 IH(MAXREF),IK(MAXREF),IL(MAXREF) INTEGER*2 IDATA(MAX3D) LOGICAL EOF, lvalms, ccpexs CHARACTER*4 KEY, CCVAL(20),KEYWORD(16) CHARACTER*400 LINE CHARACTER*80 TITLE INTEGER NTOK, IBEG(20), IEND(20), ITYP(20), IDEC(20) REAL FVAL(20) DIMENSION ADATA1(5) INTEGER IPOINT1(5) CHARACTER*30 LAB1(5),LSUSRJ(5) CHARACTER*1 CTYP1(5) REAL RSYM(4,4,64),HKL(3),CELLP(6) REAL COST(0:360),SINT(0:360) Real XPRO(MAXPRO),YPRO(MAXPRO),ZPRO(MAXPRO) REAL XR(MAXPRO),YR(MAXPRO),ZR(MAXPRO),XRT(64),YRT(64),ZRT(64) INTEGER ALPHA, BETA, GAMA REAL RRR(3,3,6),VOLL,qinv(3,3),qqq(3,3) LOGICAL ORTH COMMON /RFLX/ isolu,IP1(npeak),IP2(npeak),IP3(npeak),IP4(npeak), + IP5(npeak),IP6(npeak),FP1(npeak),FP2(npeak),FP3(npeak) DATA LAB1/'H','K','L','FP','SIGFP'/ DATA KEYWORD/'TITL','GRID','LABI','ENVI','ALPH','BETA','GAMM', + 'XRAN','YRAN','ZRAN','FLIP','RFIL','RESO','CHKO','SIGC','ORTH'/ DATA CTYP1/'H','H','H','F','Q'/ DATA IPOINT1/3*-1,2*0/ maxx = -99 maxy = -99 maxz = -99 alphamin = 0 alphamax = 180 alphastep = 3 betamin = 0 betamax = 180 betastep = 3 gammamin = 0 gammamax = 180 gammastep = 3 xmin = -99 xmax = xmin xstep = 1 ymin = -99 ymax = ymin ystep = 1 zmin = -99 zmax = zmin zstep = 1 flip = 1 rfilter = 0.6 reso = 10.0 input = 1 chkov = 0 sigcut = 0. orth = .false. ncode = 1 C ************************ CALL CCPRCS (6, 'FSEARCH', '$Date: 2002/11/04 16:34:01 $') CALL CCPFYP WRITE(6,*) WRITE(6,*)' Program: Fsearch ' WRITE(6,*)' Last update: 10/11/2001' WRITE(6,*)' Author: QUAN HAO' WRITE(6,*) DO 101 I = 0, 360 argu = 3.1415926/180.*I COST(I)=COS(argu) SINT(I)=SIN(argu) 101 CONTINUE C Set default - currently NAN CALL QNAN (VAL_MAGIC) if (ccpexs('FLMIN')) input = 1 if (ccpexs('XYZIN')) input = 2 if (ccpexs('MAPIN')) input = 3 C ******READ IN KEYWORDS ************* JJ = 0 1 LINE = ' ' NTOK = 20 CALL PARSER (KEY,LINE,IBEG,IEND,ITYP,FVAL,CCVAL, + IDEC,NTOK,EOF,.FALSE.) IF (EOF) GOTO 19 CALL CCPUPC(KEY) DO I = 1 , 16 IF (KEY.EQ.KEYWORD(I)) THEN GOTO (3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18) I END IF END DO CALL CCPERR(1,'----------- illegal keyword ------------') 3 TITLE = LINE WRITE (6,*) TITLE WRITE (6,*) GOTO 1 4 MAXX = FVAL(2) MAXY = FVAL(3) MAXZ = FVAL(4) GOTO 1 5 CALL MTZINI ITOK = 2 CALL LKYSET(LAB1,5,LSUSRJ,IPOINT1,ITOK,NTOK,LINE, + IBEG,IEND) CALL LKYIN(1,LAB1,5,NTOK,LINE,IBEG,IEND) CALL LROPEN (1,'HKLIN',1,IFAIL) IF (IFAIL.EQ.1) + CALL CCPERR(1,'------------ HKLIN failure ------------') CALL LRSYMM(1,NSYM,RSYM) CALL LRCELL(1,CELLP) CALL LRASSN (1,LAB1,5,IPOINT1,CTYP1) GOTO 1 6 input = nint(fval(2)) goto 1 7 alphamin = amod(fval(2)+360.,360.) alphamax = amod(fval(3)+360.,360.) alphastep = fval(4) goto 1 8 betamin = amod(fval(2)+360.,360.) betamax = amod(fval(3)+360.,360.) betastep = fval(4) goto 1 9 gammamin = amod(fval(2)+360.,360.) gammamax = amod(fval(3)+360.,360.) gammastep = fval(4) goto 1 10 xmin=fval(2) xmax=fval(3) xstep=fval(4) goto 1 11 ymin=fval(2) ymax=fval(3) ystep=fval(4) goto 1 12 zmin=fval(2) zmax=fval(3) zstep=fval(4) goto 1 13 flip=-1 goto 1 14 rfilter=fval(2) goto 1 15 reso=fval(2) goto 1 16 chkov=1 goto 1 17 sigcut = fval(2) goto 1 18 orth = .true. ncode = nint(fval(2)) goto 1 19 NREF = 0 MAXH = 0 MAXK = 0 MAXL = 0 C ******READ IN REFLECTION DATA************** 20 CALL LRREFF(1,RESOL,ADATA1,EOF) IF (EOF) GOTO 50 CALL IS_MAGIC (VAL_MAGIC,adata1(4),LVALMS) IF (lvalms) GOTO 20 if (adata1(4).le.0.)goto 20 IF ((1./SQRT(RESOL)).LT.reso .or. . adata1(4).lt.sigcut*adata1(5))GOTO 20 NREF = NREF + 1 DO 38 N = 1, NSYM DO 36 I = 1, 3 HKL(I) = RSYM(1,I,N)*ADATA1(1)+RSYM(2,I,N)*ADATA1(2) 1 +RSYM(3,I,N)*ADATA1(3) 36 CONTINUE MAXH = MAX(NINT(ABS(HKL(1))),MAXH) MAXK = MAX(NINT(ABS(HKL(2))),MAXK) MAXL = MAX(NINT(ABS(HKL(3))),MAXL) 38 CONTINUE IH(NREF) = NINT(ADATA1(1)) IK(NREF) = NINT(ADATA1(2)) IL(NREF) = NINT(ADATA1(3)) FF(NREF) = ADATA1(4) RES(NREF) = RESOL GOTO 20 C ******************* 50 if (xmin.eq.-99) xmin=0 if (xmax.eq.-99) xmax=cellp(1)/2 if (ymin.eq.-99) ymin=0 if (ymax.eq.-99) ymax=cellp(2)/2 if (zmin.eq.-99) zmin=0 if (zmax.eq.-99) zmax=cellp(3)/2 if (maxx.eq.-99) then do i = 1, 3 nn(i) = nint(cellp(i)) if (mod(nn(i),2).eq.1) nn(i) = nn(i) + 1 54 ntemp = nn(i) 55 ntemp1 = ntemp do j = 2,19 if (mod(ntemp,j).eq.0)ntemp = ntemp/j end do if (ntemp1.gt.ntemp) goto 55 if (ntemp.gt.19) then nn(i) = nn(i) + 2 goto 54 end if end do maxx = nn(1) maxy = nn(2) maxz = nn(3) end if IF (MAXX.LT.2*MAXH.OR.MAXY.LT.2*MAXK.OR.MAXZ.LT.2*MAXL) 1 CALL CCPERR(1, '*needs finer grids*') WRITE(6,80)MAXX,MAXY,MAXZ 80 FORMAT('Sampling intervals on xyz: ',3I8,/) IF(MAXX*MAXY*MAXZ.GT.nint(max3d*0.96)) + CALL CCPERR(1, '*out of allocated memory*') C ******************** do i = 4, 6 if(abs(cellp(i)-90.00).gt.0.01) orth = .true. end do if (orth) then call rbfro1 (cellp,voll,rrr) do i = 1,3 do j = 1,3 qinv(i,j)=rrr(i,j,ncode) end do end do call minv3(qqq,qinv,det) end if C ******************* if (input.eq.1) CALL MASK1(NUMPRO,flip,XPRO,YPRO,ZPRO) if (input.eq.2) CALL MASK2(NUMPRO,flip,XPRO,YPRO,ZPRO) if (input.eq.3) CALL MASK3(NUMPRO,flip,XPRO,YPRO,ZPRO) write(6,*)'Number of points in the mask = ',numpro write(6,*)'number of reflections = ',nref C ******MAIN LOOP BEGINS******** isolu = 0 sumr = 0. CALL CCPDPN(7,'SOLUTION','SCRATCH','U',80,0) DO 580 ALPHA = alphamin, alphamax, alphastep DO 570 BETA = betamin, betamax, betastep DO 560 GAMA = gammamin,gammamax,gammastep CALL MASK(ALPHA,BETA,GAMA,COST,SINT,XR,YR,ZR, 1 XPRO,YPRO,ZPRO,NUMPRO) C ****************** DO N = 1, NSYM C ****************** call itodata(ed,maxx,maxy,maxz,idata,n,rsym,numpro,xr,yr, 1 zr,qqq,cellp,orth,input) C *********************** call fft(ed,data,maxx,maxy,maxz,fc,nref,N, + ih,ik,il) C ********* END DO C ****Translational Shift**** DO IX = xmin,xmax,xstep DO IY = ymin,ymax,ystep DO IZ = zmin,zmax,zstep if (orth) then xrt(1)=qqq(1,1)*ix+qqq(1,2)*iy+qqq(1,3)*iz yrt(1)=qqq(2,1)*ix+qqq(2,2)*iy+qqq(2,3)*iz zrt(1)=qqq(3,1)*ix+qqq(3,2)*iy+qqq(3,3)*iz else XRT(1) = (IX)/CELLP(1) YRT(1) = (IY)/CELLP(2) ZRT(1) = (IZ)/CELLP(3) end if DO N = 2, NSYM XRT(N)=RSYM(1,1,N)*XRT(1)+RSYM(1,2,N)*YRT(1) 1 +RSYM(1,3,N)*ZRT(1) YRT(N)=RSYM(2,1,N)*XRT(1)+RSYM(2,2,N)*YRT(1) 1 +RSYM(2,3,N)*ZRT(1) ZRT(N)=RSYM(3,1,N)*XRT(1)+RSYM(3,2,N)*YRT(1) 1 +RSYM(3,3,N)*ZRT(1) END DO DO I = 1, NREF FCOMP(I)=0 DO J = 1, NSYM PH = 360*(IH(I)*XRT(J)+IK(I)*YRT(J)+IL(I)*ZRT(J)) iph = mod(nint(ph)+360000,360) SHIFT = CMPLX(COST(iPH),SINT(iPH)) FCOMP(I)=FCOMP(I)+FC(I,J)*SHIFT END DO END DO CALL RFACT(NSYM,RSYM,MAXX,MAXY,MAXZ,NREF,IH,IK,IL, 1 FF,RES,FCOMP,b,c,r) c **************************** if (chkov.eq.1) then IOVER = 0 DO JX = 0, MAXX-1 DO JY = 0, MAXY-1 DO JZ = 0, MAXZ-1 jxyz=jx*maxy*maxz+jy*maxz+jz+1 IDATA(jxyz) = 0 END DO END DO END DO DO I = 1, NUMPRO if (orth) then xrr = xr(i)+iz yrr = yr(i)+iy zrr = zr(i)+iz xrt(1)=qqq(1,1)*xrr+qqq(1,2)*yrr+qqq(1,3)*zrr yrt(1)=qqq(2,1)*xrr+qqq(2,2)*yrr+qqq(2,3)*zrr zrt(1)=qqq(3,1)*xrr+qqq(3,2)*yrr+qqq(3,3)*zrr else XRT(1) = (XR(I) + IX)/CELLP(1) YRT(1) = (YR(I) + IY)/CELLP(2) ZRT(1) = (ZR(I) + IZ)/CELLP(3) end if DO N = 2, NSYM XRT(N)=RSYM(1,1,N)*XRT(1)+RSYM(1,2,N)*YRT(1) 1 +RSYM(1,3,N)*ZRT(1)+RSYM(1,4,N) YRT(N)=RSYM(2,1,N)*XRT(1)+RSYM(2,2,N)*YRT(1) 1 +RSYM(2,3,N)*ZRT(1)+RSYM(2,4,N) ZRT(N)=RSYM(3,1,N)*XRT(1)+RSYM(3,2,N)*YRT(1) 1 +RSYM(3,3,N)*ZRT(1)+RSYM(3,4,N) END DO DO N = 1, NSYM IXRT = MOD(NINT((XRT(N)+8.)*MAXX),MAXX) IYRT = MOD(NINT((YRT(N)+8.)*MAXY),MAXY) IZRT = MOD(NINT((ZRT(N)+8.)*MAXZ),MAXZ) ixyz=ixrt*maxy*maxz+iyrt*maxz+izrt+1 IF (IDATA(ixyz).NE.N.AND. 1 IDATA(ixyz).NE.0)IOVER = IOVER + 1 IDATA(ixyz) = N END DO END DO end if C *************************************** if (r.le.rfilter) then overlap = float(iover)/(maxx*maxy*maxz) write(7)alpha, 1 beta,gama,ix,iy,iz,b,r,overlap isolu = isolu + 1 sumr = sumr + r end if END DO END DO END DO 560 CONTINUE 570 CONTINUE 580 CONTINUE if (isolu.eq.0) then write (6,*) '*** No solutions have R < ',rfilter,'***' goto 1000 end if if (chkov.eq.1) then write(6,588) 588 format (2x,'Alpha',3x,'Beta',2x,'Gamma',6x,'X',6x,'Y',6x,'Z' + ,6x,'B',8x,'R',6X,'Overlap') else write (6,589) 589 format (2x,'Alpha',3x,'Beta',2x,'Gamma',6x,'X',6x,'Y',6x,'Z' + ,6x,'B',8x,'R') end if iround = 0 600 rewind 7 if (isolu.gt.npeak) then iround = iround + 1 avr = sumr/isolu isolu = 0 sumr = 0. 650 read (7,end=600)alpha,beta,gama,ix,iy,iz,b,r,overlap if (r.gt.avr) goto 650 sumr = sumr + r isolu = isolu + 1 goto 650 else isolu = 0 660 read (7,end=680)alpha,beta,gama,ix,iy,iz,b,r,overlap if (r.gt.avr.and.iround.gt.0) goto 660 isolu = isolu + 1 ip1(isolu) = alpha ip2(isolu) = beta ip3(isolu) = gama ip4(isolu) = ix ip5(isolu) = iy ip6(isolu) = iz fp1(isolu) = b fp2(isolu) = r fp3(isolu) = overlap goto 660 680 if (isolu.gt.1) call sort(fp2,6,3) do i = 1, isolu if (chkov.eq.1)then write (6,'(6i7,2f10.3,f10.5)')IP1(I),IP2(I),IP3(I), 1 IP4(I),IP5(I),IP6(I),FP1(I),FP2(I),FP3(I) else write (6,'(6i7,2f10.3)')IP1(I),IP2(I),IP3(I), 1 IP4(I),IP5(I),IP6(I),FP1(I),FP2(I) end if end do end if 1000 CALL CCPERR(0,'Normal termination') END C----------------------------------------------------------------------- SUBROUTINE SORT(A,NI,NF) parameter(npeak=200) INTEGER NI,NF COMMON /RFLX/ NUMB,IP(npeak,6),FP(npeak,3) DIMENSION A(npeak),IP1(6),FP1(3) INT=NUMB 100 INT=INT/2 IF(2*(INT/2).EQ.INT) INT=INT-1 IFIN=NUMB-INT DO 1000 II=1,IFIN I=II J=I+INT IF(A(I).LE.A(J)) GOTO 1000 A1=A(J) DO 200 K=1,NI IP1(K)=IP(J,K) 200 CONTINUE DO 300 K=1,NF FP1(K)=FP(J,K) 300 CONTINUE 400 DO 500 K=1,NI IP(J,K)=IP(I,K) 500 CONTINUE DO 600 K=1,NF FP(J,K)=FP(I,K) 600 CONTINUE J=I I=I-INT IF(I.LE.0) GOTO 700 IF(A(I).GT.A1) GOTO 400 700 DO 800 K=1,NI IP(J,K)=IP1(K) 800 CONTINUE DO 900 K=1,NF FP(J,K)=FP1(K) 900 CONTINUE 1000 CONTINUE IF(INT.GT.1) GOTO 100 RETURN END C********************************************************************** subroutine itodata(data,maxx,maxy,maxz,idata,n,rsym,numpro, 1 xr,yr,zr,qqq,cellp,orth,input) PARAMETER(MAXPRO=200000) real data(0:maxx-1,0:maxy-1,0:maxz-1) integer*2 idata(0:maxx-1,0:maxy-1,0:maxz-1) REAL RSYM(4,4,64),qqq(3,3),cellp(6) logical orth REAL XR(MAXPRO),YR(MAXPRO),ZR(MAXPRO),XRT(64),YRT(64),ZRT(64) DO JX = 0, MAXX-1 DO JY = 0, MAXY-1 DO JZ = 0, MAXZ-1 IDATA(JX,JY,JZ) = 0 END DO END DO END DO C **************** DO I = 1, NUMPRO if (orth) then xrt(1)=qqq(1,1)*xr(i)+qqq(1,2)*yr(i)+qqq(1,3)*zr(i) yrt(1)=qqq(2,1)*xr(i)+qqq(2,2)*yr(i)+qqq(2,3)*zr(i) zrt(1)=qqq(3,1)*xr(i)+qqq(3,2)*yr(i)+qqq(3,3)*zr(i) else XRT(1) = xr(i)/CELLP(1) YRT(1) = yr(i)/CELLP(2) ZRT(1) = zr(i)/CELLP(3) end if C **************** C NB: In the case that N=1, RSYM encodes the identity operation C so the values of XRT, YRT, ZRT are unchanged. if (N.gt.1) then XRT(N)=RSYM(1,1,N)*XRT(1)+RSYM(1,2,N)*YRT(1) 1 +RSYM(1,3,N)*ZRT(1)+RSYM(1,4,N) YRT(N)=RSYM(2,1,N)*XRT(1)+RSYM(2,2,N)*YRT(1) 1 +RSYM(2,3,N)*ZRT(1)+RSYM(2,4,N) ZRT(N)=RSYM(3,1,N)*XRT(1)+RSYM(3,2,N)*YRT(1) 1 +RSYM(3,3,N)*ZRT(1)+RSYM(3,4,N) end if C ***************** IXRT = MOD(NINT((XRT(N)+8.)*MAXX),MAXX) IYRT = MOD(NINT((YRT(N)+8.)*MAXY),MAXY) IZRT = MOD(NINT((ZRT(N)+8.)*MAXZ),MAXZ) IDATA(IXRT,IYRT,IZRT) = 1 if (input.eq.2) then idata(mod(jx+1,maxx),jy,jz) = 1 idata(mod(jx-1+maxx,maxx),jy,jz) = 1 idata(jx,mod(jy+1,maxy),jz) = 1 idata(jx,mod(jy-1+maxy,maxy),jz) = 1 idata(jx,jy,mod(jz+1,maxz)) = 1 idata(jx,jy,mod(jz-1+maxz,maxz)) = 1 end if C ************** END DO DO JX = 0, MAXX-1 DO JY = 0, MAXY-1 DO JZ = 0, MAXZ-1 if (input.eq.2) then DATA(JX,JY,JZ) = idata(jx,jy,jz) else DATA(JX,JY,JZ) = idata(jx,jy,jz)+ 1 idata(mod(jx+1,maxx),jy,jz)+ 1 idata(mod(jx-1+maxx,maxx),jy,jz)+ 1 idata(jx,mod(jy+1,maxy),jz)+ 1 idata(jx,mod(jy-1+maxy,maxy),jz)+ 1 idata(jx,jy,mod(jz+1,maxz))+ 1 idata(jx,jy,mod(jz-1+maxz,maxz)) if (abs(data(jx,jy,jz)).gt.3.9) then data(jx,jy,jz)=1.0 else data(jx,jy,jz)=0.0 end if end if END DO END DO END DO return end C********************************************************************** SUBROUTINE fft(ed,data1,ix,iy,iz,fc,nref,n,ih,ik,il) real ed(0:ix-1,0:iy-1,0:iz-1) complex data1(0:ix-1,0:iy-1,0:iz-1) PARAMETER(MAXREF=120000) COMPLEX FC(MAXREF,64) INTEGER*2 IH(MAXREF),IK(MAXREF),IL(MAXREF) INTEGER h,k,l call fft3rh( ed(0,0,0),ix,iy,iz) do 10 h=0,ix-1 do 10 k=0,iy-1 do 10 l=0,iz/2-1 c *** complex part sign changed *** data1(h,k,l)=cmplx(ed(h,k,2*l),-ed(h,k,2*l+1)) if(l.ne.0) then data1(h,k,iz-l)=cmplx(ed(h,k,2*l),ed(h,k,2*l+1)) end if 10 continue DO I=1,NREF iht = mod(ih(i)+ix,ix) ikt = mod(ik(i)+iy,iy) ilt = mod(il(i)+iz,iz) FC(I,N) = DATA1(iht,ikt,ilt) c write (6,*)ih(i),ik(i),il(i),fc(i,n) END DO c write(6,*) (fc(i,n),IH(I),IK(I),IL(I), i=1,50) c write(6,*) (((data1(h,k,l),h=0,4),k=0,4),l=0,4) c stop return end C********************************************************************** c c ---------------------------------------------------------------------- c subroutine fft3rh(x,nu,nv,nw) c include 'dmheader.fh' c 3-d Real to Hermitian FFT c In: real x(u,v,w) c Out: complex x(h,k,l) (l=0...nl/2) c (complex elements are in odd l planes) c x MUST BE LARGE ENOUGH TO HOLD ((w+2)*u*v) ELEMENTS integer nu,nv,nw,nuv,nuvw real x(0:*) integer d(5) c nuv=nu*nv nuvw=nuv*(nw+2) c Real transform l to w c *** NOTE THAT THE w=nw/2 ELEMENTS ARE LOST *** d(1)=nuvw d(2)=2*nuv d(3)=nuvw d(4)=nuv d(5)=1 call realft(x(0),x(nuv),nw/2,d) c Transform k to v d(1)=nuvw d(2)=nu d(3)=2*nuv d(4)=nu d(5)=1 call cmplft(x(0),x(nuv),nv,d) c Transform h to u d(1)=nuvw d(2)=1 d(3)=2*nuv d(4)=nuv d(5)=nu call cmplft(x(0),x(nuv),nu,d) return end C********************************************************************** SUBROUTINE RFACT(NSYM,RSYM,MAXH,MAXK,MAXL,NREF,IH,IK,IL,FF, 1 RES,FCOMP,b,c,r) PARAMETER(MAXREF=120000) COMPLEX FCOMP(MAXREF) REAL RES(MAXREF),FCSUM(MAXREF),FF(MAXREF) INTEGER*2 IH(MAXREF),IK(MAXREF),IL(MAXREF) REAL RSYM(4,4,64) A1 = 0. A2 = 0. C1 = 0. C2 = 0. DO I = 1,NREF FCSUM(I)=ABS(FCOMP(I)) IF (FCSUM(I) .EQ. 0)FCSUM(I)=0.001 A1 = A1 + RES(I)**2 A2 = A2 + RES(I) TEMP = LOG(FCSUM(I)/FF(I)) c write(6,*)FCSUM(I),ih(i),ik(i),il(i),ff(i),res(i),temp C1 = C1 + TEMP*RES(I) C2 = C2 + TEMP END DO B1 = C1*A2 - C2*A1 B2 = C1*NREF - C2*A2 B3 = A1*NREF - A2*A2 B = -B2/B3 C = B1/B3 D1 = 0. D2 = 0. DO 999 I = 1, NREF FCSUM(I) = FCSUM(I)*EXP(C+B*RES(I)) D1 = D1 + ABS(FCSUM(I)-FF(I)) D2 = D2 + FF(I) 999 continue R = D1/D2 RETURN END C********************************************************************** SUBROUTINE MASK1(NUMPRO,flip,XPRO,YPRO,ZPRO) PARAMETER(MAXPRO=200000) real FLMR(0:9,0:9),FLMI(0:9,0:9),ix,iy,iz real XPRO(MAXPRO),YPRO(MAXPRO),ZPRO(MAXPRO) DATA FLMR/100*0./ DATA FLMI/100*0./ CALL ccpdpn (1,'FLMIN','OLD','F',1,IFAIL) IF (IFAIL.EQ.1) STOP lmax = 0 102 READ(1,*,END=140)L,M,F1,F2 if (l.gt.lmax) lmax=l FLMR(L,M)=F1 FLMI(L,M)=F2 GOTO 102 140 NUMPRO = 0 DO 290 IX = -60, 60 DO 280 IY = -60, 60 DO 270 IZ = -60, 60 FTP = SQRT(IX*IX+IY*IY+IZ*IZ) if(ftp.ne.0.)THETA = ACOS(IZ/FTP) if(ix.ne.0.or.iy.ne.0.)PHI = ATAN2(IY,IX) X = COS(THETA) FTPR = 0.0 DO 260 L = 0, lmax DO 250 M = 0, L IF(ABS(FLMR(L,M)).LT.0.00001.AND. * ABS(FLMI(L,M)).LT.0.00001)GOTO 250 PLM = PLGNDR(L,M,X) FTPR = FTPR + FLMR(L,M)*PLM*COS(M*PHI)+ * FLMI(L,M)*PLM*SIN(M*PHI) 250 CONTINUE 260 CONTINUE IF (FTP.LT.FTPR) THEN NUMPRO = NUMPRO + 1 XPRO(NUMPRO) = IX c *** model flipped over if flip = -1 *** YPRO(NUMPRO) = flip*IY ZPRO(NUMPRO) = flip*IZ END IF 270 CONTINUE 280 CONTINUE 290 CONTINUE END C C********************************************************************** SUBROUTINE MASK2(NUMPRO,flip,XPRO,YPRO,ZPRO) PARAMETER(MAXPRO=200000) character*3 lett real XPRO(MAXPRO),YPRO(MAXPRO),ZPRO(MAXPRO) CALL ccpdpn (2,'XYZIN','OLD','F',1,IFAIL) IF (IFAIL.EQ.1) STOP 1 read (2,'(A3)') lett if (lett.ne.'ATO') goto 1 backspace 2 NUMPRO = 0 10 read (2,'(A3,27X,3F8.3)')lett,x,y,z if (lett.eq.'END'.or.lett.eq.' '.or.lett.eq.'EOF') goto 100 if (lett.ne.'ATO') goto 10 NUMPRO = NUMPRO + 1 XPRO(NUMPRO) = X c *** model flipped over if flip = -1 *** YPRO(NUMPRO) = flip*Y ZPRO(NUMPRO) = flip*Z goto 10 100 continue END C C********************************************************************** SUBROUTINE MASK3(NUMPRO,flip,XPRO,YPRO,ZPRO) PARAMETER(MAXPRO=200000) PARAMETER (NUVLIM=160000) real XPRO(MAXPRO),YPRO(MAXPRO),ZPRO(MAXPRO) REAL RHMAX,RHMEAN,RHMIN,RHRMS INTEGER LMODE,LSPGRP,NSEC,NU1,NU2,NV1,NV2,NW1 CHARACTER*80 TITLE REAL CELL(6),XYZPRO(3) INTEGER IUVW(3),MXYZ(3) DIMENSION RHOSEC(NUVLIM),rrr(3,3,6) ncode = 1 CALL MRDHDR(3,'MAPIN',TITLE,NSEC,IUVW,MXYZ,NW1,NU1,NU2, + NV1,NV2,CELL,LSPGRP,LMODE,RHMIN,RHMAX,RHMEAN, + RHRMS) if((nu2-nu1+1)*(nv2-nv1+1).gt.nuvlim)stop + '*** ARRAY SIZE NOT BIG ENOUGH ***' call rbfro1 (cell,voll,rrr) NUMPRO = 0 do lsec=nw1,nw1+nsec-1 call mposn(3,lsec) call mgulpr(3,rhosec,ier) if(ier.ne.0)stop '*** Map Reading Error ***' do j=1,nv2-nv1+1 do i=1,nu2-nu1+1 if(rhosec((j-1)*(nv2-nv1+1)+i).gt.0.1*RHRMS) then NUMPRO = NUMPRO + 1 ip = iuvw(1) xyzpro(ip)=float(i-1)/mxyz(ip) ip = iuvw(2) xyzpro(ip)=float(j-1)/mxyz(ip) ip = iuvw(3) xyzpro(ip)=float(lsec)/mxyz(ip) XPRO(NUMPRO) = rrr(1,1,ncode)*xyzpro(1) + +rrr(1,2,ncode)*xyzpro(2) + rrr(1,3,ncode)*xyzpro(3) c *** model flipped over if flip = -1 *** YPRO(NUMPRO) = flip*(rrr(2,1,ncode)*xyzpro(1) + +rrr(2,2,ncode)*xyzpro(2) + rrr(2,3,ncode)*xyzpro(3)) ZPRO(NUMPRO) = flip*(rrr(3,1,ncode)*xyzpro(1) + +rrr(3,2,ncode)*xyzpro(2) + rrr(3,3,ncode)*xyzpro(3)) end if end do end do end do END C C ****************************** C Calculate Legendre polynomials C ****************************** FUNCTION plgndr(l,m,x) INTEGER l,m REAL plgndr,x INTEGER i,ll REAL fact,pll,pmm,pmmp1,somx2 if ((m.lt.0).or.(m.gt.l).or.(abs(x).gt.1.)) then write(6,*) 'bad arguments in plgndr' stop end if pmm=1. if(m.gt.0) then somx2=sqrt((1.-x)*(1.+x)) fact=1. do 11 i=1,m pmm=-pmm*fact*somx2 fact=fact+2. 11 continue endif if(l.eq.m) then plgndr=pmm else pmmp1=x*(2*m+1)*pmm if(l.eq.m+1) then plgndr=pmmp1 else do 12 ll=m+2,l pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m) pmm=pmmp1 pmmp1=pll 12 continue plgndr=pll endif endif return END C *********************************** C Rotate mask C *********************************** SUBROUTINE MASK(ALPHA,BETA,GAMA,COST,SINT,XR,YR,ZR, 1 XPRO,YPRO,ZPRO,NUMPRO) PARAMETER(MAXPRO=200000) real XPRO(MAXPRO),YPRO(MAXPRO),ZPRO(MAXPRO) REAL XR(MAXPRO),YR(MAXPRO),ZR(MAXPRO),ix, iy, iz REAL COST(0:360),SINT(0:360) INTEGER ALPHA, BETA, GAMA DO 290 I = 1, NUMPRO IX = XPRO(I) IY = YPRO(I) IZ = ZPRO(I) XR(I) = (COST(ALPHA)*COST(BETA)*COST(GAMA)- * SINT(ALPHA)*SINT(GAMA))*IX + * (-COST(ALPHA)*COST(BETA)*SINT(GAMA)- * SINT(ALPHA)*COST(GAMA))*IY + * COST(ALPHA)*SINT(BETA)*IZ YR(I) = (SINT(ALPHA)*COST(BETA)*COST(GAMA)+ * COST(ALPHA)*SINT(GAMA))*IX + * (-SINT(ALPHA)*COST(BETA)*SINT(GAMA)+ * COST(ALPHA)*COST(GAMA))*IY + * SINT(ALPHA)*SINT(BETA)*IZ ZR(I) = -SINT(BETA)*COST(GAMA)*IX + * SINT(BETA)*SINT(GAMA)*IY + COST(BETA)*IZ 290 CONTINUE END