C ++******************************************************************** C * C POLQS.F * C OPFILEC FEB 03 ARDEAN LEITH * C FILELIST FEB 06 ARDEAN LEITH * C * C ********************************************************************** C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM. AUTHOR: J.FRANK * C=* Copyright (C) 1985-2007 Health Research Inc. * C=* * C=* HEALTH RESEARCH INCORPORATED (HRI), * C=* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455. * C=* * C=* Email: spider@wadsworth.org * C=* * C=* This program is free software; you can redistribute it and/or * C=* modify it under the terms of the GNU General Public License as * C=* published by the Free Software Foundation; either version 2 of the * C=* License, or (at your option) any later version. * C=* * C=* This program is distributed in the hope that it will be useful, * C=* but WITHOUT ANY WARRANTY; without even the implied warranty of * C=* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * C=* General Public License for more details. * C=* * C=* You should have received a copy of the GNU General Public License * C=* along with this program; if not, write to the * C=* Free Software Foundation, Inc., * C=* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. * C=* * C ********************************************************************** C * C POLQS(MAXMEM) * C * C PURPOSE: * C * C PARAMETERS: * C * C 02/22/07 USED FILELIST C C 07/14/05 Corrected minor bug with 1D projection C OpenMP switched off - it does not work on Linux C C 07/29/99 PARALLEL, F90 VERSION; CORRECTED ON 02/13/02 C C 02/11/99 TRIGONOMETRIC FUNCTIONS CHANGED TO F90. C C 08/28/97 C DISTRIBUTION OF ANGLES IN ANGTAB CHANGED TO AGREE WITH CURRENT VO EA C 10/3/96 C COMPLETE VERSION FOR SPIDER C CALCULATES THE 1D PROJECTIONS OF 2D IMAGES (NOT NORMALIZED). C 11/19/90 C 05/16/95 C SAME AS POLQT, BUT FOURIER VERSION ... C WEIGHTS ON CIRCLE AS IN POLQW C PROBLEM OF DUPLICATED LINES IN WEIGHTING SOLVED C WEIGHTS TAKEN OUT OF THE DO-LOOP C CORRECT LINE CALCULATION, EXCLUSION OF IDENTICAL DIRECTIONS ... C C ********************************************************************** SUBROUTINE POLQS(MAXMEM) PARAMETER (NILMAX=4000) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' INCLUDE 'F90ALLOC.INC' CHARACTER*80 FINPAT,FINPIC,FINFO,FINDOC COMMON /F_SPEC/ FINPAT,FINPIC,FINFO,FINDOC,NLET CHARACTER*1 MODIS COMMON /POLS/ MODIS COMMON BUF(1024) COMMON /PAR/ LDP,NM,LDPNM DIMENSION TARRAY(2) REAL, ALLOCATABLE, DIMENSION(:,:,:) :: PRJ REAL, ALLOCATABLE, DIMENSION(:,:) :: P REAL, DIMENSION(:,:), POINTER :: IPQ DATA INPIC/99/,INDOC/96/,NDOC/95/ C N - LINEAR DIMENSION OF 2D IMAGE C NANG - NUMBER OF ANGLES C NIMA - NUMBER OF 2D IMAGES NMAX = NIMAX CALL FILELIST(.TRUE.,NDOC,FINPAT,NLET,INUMBR,NMAX,NIMA, & 'INPUT FILE TEMPLATE (E.G. PIC****)',IRTFLG) IF (IRTFLG .NE. 0) RETURN C GET FIRST PICTURE TO DETERMINE DIMS CALL FILGET(FINPAT,FINPIC,NLET,INUMBR(1),IRTFLG) IF (IRTFLG .NE. 0) RETURN MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINPIC,INPIC,'O',IFORM, & N,NROW,NSL,MAXIM,' ',.FALSE.,IRTFLG) IF (IRTFLG.NE.0) RETURN CLOSE(INPIC) C NANG - TOTAL NUMBER OF ANGLES CALL RDPRMI(LENPROJ,NANG,NOT_USED,'LENGTH OF LINE PROJECTION') CALL RDPRM2(F1,F2,NOT_USED,'MINIMUM, MAXIMUM FREQUENCY') IF (F1.LT.0.0 .OR. F2.GT.0.5 .OR. F1.GT.F2) THEN CALL ERRT(14,'OP',IER) RETURN ENDIF IF (F1.EQ.0.0 .AND. F2.EQ.0.0) THEN F1=0.0 F2=0.5 ENDIF CALL RDPRMI(IDPSI,NANG,NOT_USED, 'ACCURACY OF THETA') C PAP 07/14/05 C I cannot fix the code for arbitrary NANG (number of psi angles). C The reason is that in too many places the code depends on the assumption that NANG=180. C The main problem is in SPIN, where the spin step is IDPSI, which again can be any... C Since all angles are coded by their integer numbers (i.e., l1=3 means angle 3*180/NANG) C the corrections are virtually impossible without major rewriting of the code. C IF (NANG.EQ.0) THEN C DEFAULT 1 DEGREE ACCURACY FOR PSI NANG=180 C ELSE C IDUM=180/NANG C IF(IDUM*NANG.NE.180) THEN C WRITE(NOUT,*) ' 180 has to be divisible by NANG' C CALL ERRT(14,'OP',NE) C GOTO 991 C ENDIF C NANG=IDUM C ENDIF C DEFAULT 5 DEGREES ACCURACY FOR THETA IF (IDPSI.EQ.0) IDPSI=5 CALL RDPRMI(MAXIT,IDUM,NOT_USED, 'MAXIMUM NUMBER OF CYCLES') DPSI=IDPSI CALL GENAT(BUF,DPSI,NANGT,.FALSE.) MODIS = 'E' RI = LENPROJ/2 LENPROJ = LENPROJ/2 LENPROJ = LENPROJ*2+1 IF (LENPROJ .GT. N-2) THEN LENRPOJ = N-2 IF (MOD(N,2) .EQ. 0) LENRPOJ = LENRPOJ-1 WRITE(NOUT,*) & 'LENGTH OF PROJECTION TOO LARGE! LIMITED TO ',LENRPOJ ENDIF C ALWAYS SKIP (0) TERM, LENB has to be odd to fall on the real part LENB=MAX(NINT(2*F1*LENPROJ+1),3) LENB=LENB+MOD(LENB+1,2) LENE=MIN(NINT(2*F2*LENPROJ),LENPROJ-1) LENF=LENE-LENB+1 LDP=N/2+1 LDPNM=LENPROJ/2+1 ALLOCATE(PRJ(LENF,2*NANG,NIMA),STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'OP, PRJ',IER) RETURN ENDIF TM1 = ETIME(TARRAY) CALL RDL_P(N,NANG,LENPROJ,RI,INUMBR,NIMA,PRJ,LENB,LENF) TM2 = ETIME(TARRAY) WRITE(NOUT,*) ' TIME FOR READING PROJECTIONS = ',TM2-TM1 MAXKEY = NIMA MAXREG = 4 CALL GETDOCDAT('DOCUMENT WITH INITIAL ANGLES',.TRUE.,FINDOC, & INDOC,.FALSE.,MAXREG,MAXKEY,IPQ,IRTFLG) IF (IRTFLG.NE.0) RETURN ALLOCATE(P(MAXREG-1,NIMA),STAT=IRTFLG) IF (IRTFLG.NE.0) CALL ERRT(46,'OP, P',IER) DO I=1,MAXKEY DO J=1,3 P(J,I) = IPQ(J+1,I) ENDDO ENDDO C IPQ NO LONGER NEEDED DEALLOCATE(IPQ) TM3 = ETIME(TARRAY) CALL FCANG(NANG,NIMA,LENF,MAXIT,P,PRJ,NANGT,IDPSI) TM4 = ETIME(TARRAY) WRITE(NOUT,*) ' TIME FOR ANGLE SEARCH = ',TM4-TM3 991 IF(ALLOCATED(P)) DEALLOCATE(P) IF(ALLOCATED(PRJ)) DEALLOCATE(PRJ) END SUBROUTINE FCANG(NANG,NIMA,LENPROJ,MAXIT, & P,PRJ,NANGT,IDPSI) INCLUDE 'CMBLOCK.INC' PARAMETER (NDLI=4) DIMENSION DLIST(NDLI) DIMENSION P(3,NIMA) DIMENSION PRJ(LENPROJ,2*NANG,NIMA) REAL, ALLOCATABLE, DIMENSION(:,:) :: ANGTAB CHARACTER*8 CTIME DOUBLE PRECISION DIS(NIMA) DOUBLE PRECISION DIST DATA NDOC/77/ C NEW TABLES. C ANGTAB(1 - THETA C ANGTAB(2 - PHI LOGICAL CHANGE,CHANGEPOS C SET FIRST THREE ANGLES TO ZERO FOREVER ... C P(1,1)=0.0 C P(2,1)=0.0 C P(3,1)=0.0 ALLOCATE(ANGTAB(2,NANGT),STAT=IRTFLG) IF (IRTFLG.NE.0) CALL ERRT(46,'OP, ANGTAB',IER) DPSI=IDPSI CALL GENAT(ANGTAB,DPSI,NANGT,.TRUE.) FVALUE=DIST(NIMA,NANG,LENPROJ,PRJ,P) WRITE(NOUT,2177) FVALUE 2177 FORMAT(' INITIAL DISCREPANCY ',1PE10.3) WRITE(NOUT,2178) 2178 FORMAT(' INITIAL ANGLES',5X,'PSI',7X,'THETA',5X,'PHI') C TRY TO SPIN ITER=0 507 CHANGE=.TRUE. CHANGEPOS=.FALSE. ITER=ITER+1 DO I=1,NIMA WRITE(NOUT,2179) I,(P(J,I),J=1,3) ENDDO 2179 FORMAT(I15,3(3X,F7.2)) CALL MYTIME(CTIME) WRITE(NOUT,777) CTIME 777 FORMAT(' ENTERED SPIN : ',A8) DO ISPIN=1,NIMA C CALL SPINSPIN C & (NIMA,NANG,LENPROJ,PRJ,P, C & ISPIN,IDPSI,ANGTAB,NANGT,IPSIM,IANGMIN,CHANGE, C & L1,L2,LT,LT1,LT2,LORD,W) CALL SPIN(NIMA,NANG,LENPROJ,PRJ,P, & ISPIN,IDPSI,ANGTAB,NANGT,IPSIM,IANGMIN,CHANGE) IF (CHANGE) THEN CHANGEPOS=.TRUE. P(1,ISPIN)=IPSIM P(2,ISPIN)=ANGTAB(1,IANGMIN) P(3,ISPIN)=ANGTAB(2,IANGMIN) ENDIF ENDDO IF (CHANGEPOS) THEN FVALUE=DIST(NIMA,NANG,LENPROJ,PRJ,P) WRITE(NOUT,2181) ITER,FVALUE 2181 FORMAT(' CYCLE',I9,' NEW VALUE=',1PE10.3) DLIST(1)=-1 DLIST(2)=ITER DLIST(3)=FVALUE CALL SAVD(NDOC,DLIST,3,IRTFLG) DO I=1,NIMA DLIST(1)=I DO J=1,3 DLIST(J+1)=P(J,I) ENDDO CALL SAVD(NDOC,DLIST,NDLI,IRTFLG) ENDDO ENDIF IF (CHANGEPOS.AND.(MAXIT.EQ.0.OR.ITER.LT.MAXIT)) GOTO 507 CLOSE(NDOC) CALL SAVDC DEALLOCATE(ANGTAB) END SUBROUTINE SPIN(NIMA,NANG,LENPROJ,PRJ,P, & ISPIN,IDPSI,ANGTAB,NANGT,MMIN,LMIN,CHANGE) DIMENSION P(3,NIMA),PRJ(LENPROJ,2*NANG,NIMA) DIMENSION ANGTAB(2,NANGT) DOUBLE PRECISION DMI(NANGT),DMIN DIMENSION LMI(NANGT),MMI(NANGT) C,EDL1 LOGICAL CHANGE C LENGTH OF ARRAYS NEEDED IN DWTG LNM2=NIMA*(NIMA-1)/2 C SPIN Cc$omp parallel do private(l) DO L=1,NANGT CALL DWTG(P,ANGTAB(1,L),PRJ,ISPIN,NIMA,NANG,LENPROJ, & L,IDPSI,DMI(L),LMI(L),MMI(L),LNM2) ENDDO DMIN=1.0D23 LMIN=-1 MMIN=-1 DO 2 L=1,NANGT DO I=1,NIMA IF(I.NE.ISPIN.AND. & P(2,I).EQ.ANGTAB(1,L).AND.P(3,I).EQ.ANGTAB(2,L)) GOTO 2 ENDDO IF (DMI(L).LT.DMIN) THEN DMIN=DMI(L) LMIN=LMI(L) MMIN=MMI(L) ENDIF 2 CONTINUE C CHECK WHETHER POSITION CHANGED C PRINT *,NANGT,DMIN,LMIN,MMIN C PRINT *,FLOAT(MMIN),ANGTAB(1,LMIN),ANGTAB(2,LMIN) IF (P(1,ISPIN).EQ.FLOAT(MMIN)) THEN IF (P(2,ISPIN).EQ.ANGTAB(1,LMIN)) THEN IF (P(3,ISPIN).EQ.ANGTAB(2,LMIN)) THEN CHANGE = .FALSE. RETURN ENDIF ENDIF ENDIF CHANGE = .TRUE. END SUBROUTINE DWTG(PP,ANGTAB,PRJ,ISPIN,NIMA,NANG,LENPROJ, & L,IDPSI,DMIN,LMIN,MMIN,LNM2) DIMENSION PP(3,NIMA),ANGTAB(2),PRJ(LENPROJ,2*NANG,NIMA) C AUTOMATIC ARRAYS DIMENSION L1(LNM2),L2(LNM2),LT(LNM2) DIMENSION W(LNM2),LT1(LNM2),LT2(LNM2),LORD(LNM2) DIMENSION P(3,NIMA) DOUBLE PRECISION DMIN,DID,DIST DMIN=1.0D23 LMIN=-1 MMIN=-1 DO I=1,NIMA IF(I.NE.ISPIN.AND. & PP(2,I).EQ.ANGTAB(1).AND.PP(3,I).EQ.ANGTAB(2)) RETURN ENDDO P=PP NIM=NIMA-1 P(2,ISPIN)=ANGTAB(1) P(3,ISPIN)=ANGTAB(2) P(1,ISPIN)=0.0 K=0 DO I=1,NIMA-1 DO J=I+1,NIMA IF(I.EQ.ISPIN.OR.J.EQ.ISPIN) THEN K=K+1 LORD(K)=K CALL CALDILW1(P(1,I),P(1,J),NANG,L1(K),L2(K)) ENDIF ENDDO ENDDO K=0 DO I=1,NIMA-1 DO J=I+1,NIMA IF(I.EQ.ISPIN) THEN K=K+1 LT(K)=L1(K) IF(LT(K).GT.NANG) LT(K)=LT(K)-NANG ELSEIF(J.EQ.ISPIN) THEN K=K+1 LT(K)=L2(K) IF(LT(K).GT.NANG) LT(K)=LT(K)-NANG ENDIF ENDDO ENDDO C FROM SMALL TO LARGE CALL SORT2(LT,LORD,NIM) QT=1.0 DO K=2,NIM-1 IF (LT(K).NE.LT(K-1)) THEN W(LORD(K))=FLOAT(LT(K+1)-LT(K-1))/2 QT=1.0 ELSE QT=QT+1.0 W(LORD(K))=-QT ENDIF ENDDO IF (LT(NIM).NE.LT(NIM-1)) THEN W(LORD(NIM))=FLOAT(LT(1)+NANG-LT(NIM-1))/2 ELSE QT=QT+1 W(LORD(NIM))=-QT ENDIF W(LORD(1))=FLOAT(LT(2)-LT(NIM)+NANG)/2 K=1 105 K=K+1 IF (W(LORD(K)).LT.0.0) THEN KB=K-1 K=KB 106 K=K+1 IF(W(LORD(K)).LT.0.0) THEN IF(K.NE.NIM) GOTO 106 ELSE K=K-1 ENDIF W(LORD(KB))=-W(LORD(KB))/W(LORD(K)) DO KK=KB+1,K W(LORD(KK))=W(LORD(KB)) ENDDO ENDIF IF(K.LT.NIM) GOTO 105 WT=0.0 DO K=1,NIM WT=WT+W(K) ENDDO DO K=1,NIM W(K)=W(K)/WT ENDDO DO M=0,359,IDPSI P(1,ISPIN) = M C GET LINES ROTATED BY PSI ANGLE K=0 DO I=1,NIMA-1 DO J=I+1,NIMA IF (I.EQ.ISPIN.OR.J.EQ.ISPIN) THEN K=K+1 IF (L1(K).EQ.NANG/2+1.OR.L1(K).EQ.271 .OR. & L2(K).EQ.91.OR.L2(K).EQ.271) THEN CALL CALDILW1(P(1,I),P(1,J),NANG,LT1(K),LT2(K)) ELSE IF(I.EQ.ISPIN) THEN LT1(K)=MOD(L1(K)-M+2*NANG-1,2*NANG)+1 LT2(K)=L2(K) ELSE LT1(K)=L1(K) LT2(K)=MOD(L2(K)-M+2*NANG-1,2*NANG)+1 ENDIF ENDIF ENDIF ENDDO ENDDO C DO THE DISTANCE WITH WEIGHTING DIST=0.0D0 K=0 DO I=1,NIMA-1 DO J=I+1,NIMA IF (I.EQ.ISPIN.OR.J.EQ.ISPIN) THEN K=K+1 if (w(k).le.0.0) then print *,' incorrect weight',k,w(k) stop endif DIST = DIST+EDL1(PRJ(1,LT1(K),I), & PRJ(1,LT2(K),J),LENPROJ)*W(K) ENDIF ENDDO ENDDO IF (DIST.LT.DMIN) THEN DMIN=DIST LMIN=L MMIN=M ENDIF ENDDO C END OF DO-LOOP OVER M (PSI ANGLE) END SUBROUTINE SPINSPIN(NIMA,NANG,LENPROJ,PRJ,P, & ISPIN,IDPSI,ANGTAB,NANGT,MMIN,LMIN,CHANGE, & L1,L2,LT,LT1,LT2,LORD,W) DIMENSION P(3,NIMA),PRJ(LENPROJ,2*NANG,NIMA),PTMP(3) DIMENSION ANGTAB(2,NANGT) DOUBLE PRECISION DMIN,DID,DIST C,EDL1 LOGICAL CHANGE DIMENSION L1(*),L2(*),LT(*),LT1(*),LT2(*),LORD(*) DIMENSION W(*) NIM=NIMA-1 C SPIN DMIN=1.0E23 LMIN=-1 MMIN=-1 DO 1 L=1,NANGT DO I=1,NIMA IF(I.NE.ISPIN.AND. & P(2,I).EQ.ANGTAB(1,L).AND.P(3,I).EQ.ANGTAB(2,L)) GOTO 1 ENDDO PTMP(2)=P(2,ISPIN) PTMP(3)=P(3,ISPIN) PTMP(1)=P(1,ISPIN) P(2,ISPIN)=ANGTAB(1,L) P(3,ISPIN)=ANGTAB(2,L) P(1,ISPIN)=0.0 K=0 DO 14 I=1,NIMA-1 DO 14 J=I+1,NIMA IF (I.EQ.ISPIN.OR.J.EQ.ISPIN) THEN K=K+1 LORD(K)=K CALL CALDILW1(P(1,I),P(1,J),NANG,L1(K),L2(K)) ENDIF 14 CONTINUE K=0 DO 101 I=1,NIMA-1 DO 101 J=I+1,NIMA IF (I.EQ.ISPIN) THEN K=K+1 LT(K)=L1(K) IF (LT(K).GT.NANG) LT(K)=LT(K)-NANG ELSEI F(J.EQ.ISPIN) THEN K=K+1 LT(K)=L2(K) IF(LT(K).GT.NANG) LT(K)=LT(K)-NANG ENDIF 101 CONTINUE C FROM SMALL TO LARGE CALL SORT2(LT,LORD,NIM) QT=1.0 DO 102 K=2,NIM-1 IF(LT(K).NE.LT(K-1)) THEN W(LORD(K))=FLOAT(LT(K+1)-LT(K-1))/2 QT=1.0 ELSE QT=QT+1.0 W(LORD(K))=-QT ENDIF 102 CONTINUE IF(LT(NIM).NE.LT(NIM-1)) THEN W(LORD(NIM))=FLOAT(LT(1)+NANG-LT(NIM-1))/2 ELSE QT=QT+1 W(LORD(NIM))=-QT ENDIF W(LORD(1))=FLOAT(LT(2)-LT(NIM)+NANG)/2 C K=1 105 K=K+1 IF (W(LORD(K)).LT.0.0) THEN KB=K-1 K=KB 106 K=K+1 IF (W(LORD(K)).LT.0.0) THEN IF(K.NE.NIM) GOTO 106 ELSE K=K-1 ENDIF W(LORD(KB))=-W(LORD(KB))/W(LORD(K)) DO 107 KK=KB+1,K 107 W(LORD(KK))=W(LORD(KB)) ENDIF IF (K.LT.NIM) GOTO 105 WT=0.0 DO K=1,NIM WT=WT+W(K) ENDDO DO K=1,NIM W(K)=W(K)/WT ENDDO DO 2 M=0,359,IDPSI P(1,ISPIN)=M C GET LINES ROTATED BY PSI ANGLE K=0 DO I=1,NIMA-1 DO J=I+1,NIMA IF(I.EQ.ISPIN.OR.J.EQ.ISPIN) THEN K=K+1 IF(L1(K).EQ.91.OR.L1(K).EQ.271.OR.L2(K).EQ.91.OR.L2(K).EQ.271) & THEN CALL CALDILW1(P(1,I),P(1,J),NANG,LT1(K),LT2(K)) ELSE IF(I.EQ.ISPIN) THEN LT1(K)=MOD(L1(K)-M+2*NANG-1,2*NANG)+1 LT2(K)=L2(K) ELSE LT1(K)=L1(K) LT2(K)=MOD(L2(K)-M+2*NANG-1,2*NANG)+1 ENDIF ENDIF ENDIF ENDDO ENDDO C DO THE DISTANCE WITH WEIGHTING DIST=0.0D0 K=0 DO I=1,NIMA-1 DO J=I+1,NIMA IF(I.EQ.ISPIN.OR.J.EQ.ISPIN) THEN K=K+1 DID=EDL1(PRJ(1,LT1(K),I),PRJ(1,LT2(K),J),LENPROJ) DIST=DIST+DID*W(K) ENDIF ENDDO ENDDO C IF(DIST.LT.DMIN) THEN DMIN=DIST LMIN=L MMIN=M ENDIF 2 CONTINUE P(2,ISPIN)=PTMP(2) P(3,ISPIN)=PTMP(3) P(1,ISPIN)=PTMP(1) 1 CONTINUE C CHECK WHETHER POSITION CHANGED C PRINT *,NANGT,DMIN,LMIN,MMIN C PRINT *,FLOAT(MMIN),ANGTAB(1,LMIN),ANGTAB(2,LMIN) IF(P(1,ISPIN).EQ.FLOAT(MMIN)) THEN IF(P(2,ISPIN).EQ.ANGTAB(1,LMIN)) THEN IF(P(3,ISPIN).EQ.ANGTAB(2,LMIN)) THEN CHANGE=.FALSE. RETURN ENDIF ENDIF ENDIF CHANGE=.TRUE. END DOUBLE PRECISION FUNCTION DIST & (NIMA,NANG,LENPROJ,PRJ,P) DIMENSION P(3,NIMA),PRJ(LENPROJ,2*NANG,NIMA) DOUBLE PRECISION DIN,DIN1,DIN2,DID DIN=0.0D0 C#ifdef SP_MP C NTR=NIMA*(NIMA-1)/2 Cc$omp parallel do private(k,i,j,did) Cc$omp+reduction(+:din) C DO K=1,NTR C CALL DECO(K,NIMA,I,J) C CALL CALDIL C & (PRJ(1,1,I),P(1,I),PRJ(1,1,J),P(1,J),LENPROJ,NANG,DID) C DIN=DIN+DID C#else DO I=1,NIMA-1 DO J=I+1,NIMA CALL CALDIL & (PRJ(1,1,I),P(1,I),PRJ(1,1,J),P(1,J),LENPROJ,NANG,DID) DIN=DIN+DID ENDDO C#endif ENDDO DIST=DIN END SUBROUTINE DECO(M,N,I,J) CCCCCCCCCCCCCCCCCCCCCCCCCCC IFI2 !! C SUCH CHANGE MEANS THAT ANGLES FOUND SHOULD BE REVERSED TOO: C ALPHA1<->ALPHA2 DIMENSION PRJ1(LENPROJ,2*NANG),FI1(3), & PRJ2(LENPROJ,2*NANG),FI2(3) C DOUBLE PRECISION EDL1,EDL2,EDL3,DID,R1(3,3),R2(3,3),R3(3,3) DOUBLE PRECISION EDL2,EDL3,DID,R1(3,3),R2(3,3),R3(3,3) DOUBLE PRECISION ECL1,ECL2,ECL3 DOUBLE PRECISION PSI,THETA,PHI,DEPS CHARACTER*1 MODIS COMMON /POLS/ MODIS DOUBLE PRECISION QUADPI,RAD_TO_DGR PARAMETER (QUADPI = 3.141592653589793238462643383279502884197) PARAMETER (RAD_TO_DGR = (180.0/QUADPI)) PARAMETER (DEPS = 1.0D-7) CALL BLDR(R1,-FI1(3),-FI1(2),-FI1(1)) CALL BLDR(R2,FI2(1),FI2(2),FI2(3)) DO 1 I=1,3 DO 1 J=1,3 R3(I,J)=0.0 DO 1 K=1,3 1 R3(I,J)=R3(I,J)+R2(I,K)*R1(K,J) C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C LIMIT PRECISION DO 5 J=1,3 DO 5 I=1,3 IF(DABS(R3(I,J)).LT.DEPS) R3(I,J)=0.0D0 IF(R3(I,J)-1.0D0.GT.-DEPS) R3(I,J)=1.0D0 IF(R3(I,J)+1.0D0.LT.DEPS) R3(I,J)=-1.0D0 5 CONTINUE IF (R3(3,3).EQ.1.0) THEN THETA=0.0 PSI=0.0 IF (R3(1,1).EQ.0.0) THEN PHI=RAD_TO_DGR*DASIN(R3(1,2)) ELSE PHI=RAD_TO_DGR*DATAN2(R3(1,2),R3(1,1)) ENDIF ELSEIF(R3(3,3).EQ.-1.0) THEN THETA=180.0 PSI=0.0 IF (R3(1,1).EQ.0.0) THEN PHI=RAD_TO_DGR*DASIN(-R3(1,2)) ELSE PHI=RAD_TO_DGR*DATAN2(-R3(1,2),-R3(1,1)) ENDIF ELSE THETA=RAD_TO_DGR*DACOS(R3(3,3)) ST=DSIGN(1.0D0,THETA) IF (R3(3,1).EQ.0.0) THEN IF(ST.NE.DSIGN(1.0D0,R3(3,2))) THEN PHI=270.0 ELSE PHI=90.0 ENDIF ELSE PHI=RAD_TO_DGR*DATAN2(R3(3,2)*ST,R3(3,1)*ST) ENDIF IF (R3(1,3).EQ.0.0) THEN IF(ST.NE.DSIGN(1.0D0,R3(2,3))) THEN PSI=270.0 ELSE PSI=90.0 ENDIF ELSE PSI=RAD_TO_DGR*DATAN2(R3(2,3)*ST,-R3(1,3)*ST) ENDIF ENDIF IF (PSI.LT.0.0) PSI=PSI+360.0 IF (THETA.LT.0.0) THETA=THETA+360.0 IF (PHI.LT.0.0) PHI=PHI+360.0 C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C IF(DABS(R3(3,1)).LT.DEPS .AND. DABS(R3(3,2)).LT.DEPS) THEN C PRINT *,' ALPHA 1 CANNOT BE CALCULATED' C DID=2.0 C RETURN C ENDIF C IF(DABS(R3(1,3)).LT.DEPS .AND. DABS(R3(2,3)).LT.DEPS) THEN C PRINT *,' ALPHA 2 CANNOT BE CALCULATED' C DID=2.0 C RETURN C ENDIF ALPHA1=90.0+PHI ALPHA2=90.0-PSI IF (ALPHA1.GE.0.0) THEN L1=MOD(INT(NANG*ALPHA1/180.0),2*NANG)+1 ELSE L1=MOD(INT(NANG*(360.0+ALPHA1)/180.0),2*NANG)+1 ENDIF IF (ALPHA2.GE.0.0) THEN L2=MOD(INT(NANG*ALPHA2/180.0),2*NANG)+1 ELSE L2=MOD(INT(NANG*(360.0+ALPHA2)/180.0),2*NANG)+1 ENDIF IF (MODIS.EQ.'E') THEN DID = EDL1(PRJ1(1,L1),PRJ2(1,L2),LENPROJ) ELSEIF (MODIS.EQ.'C') THEN IF (L1.GT.NANG) THEN L1=L1-NANG IF(L2.GT.NANG) THEN L2=L2-NANG DID=ECL1(PRJ1(1,L1),PRJ2(1,L2),LENPROJ) ELSE DID=ECL2(PRJ1(1,L1),PRJ2(1,L2),LENPROJ) ENDIF ELSE IF (L2.GT.NANG) THEN L2=L2-NANG DID=ECL3(PRJ1(1,L1),PRJ2(1,L2),LENPROJ) ELSE DID=ECL1(PRJ1(1,L1),PRJ2(1,L2),LENPROJ) ENDIF ENDIF ELSE STOP ENDIF C PRINT *,ALPHA1,L1,ALPHA2,L2,DID END SUBROUTINE CALDILW1(FI1,FI2,NANG,L1,L2) C CALCULATES DISTANCE DID USING COMMON LINE C WARNING - THE RESULT IS NOT THE SAME AFTER CHANGING THE ORDER C OF ARGUMENTS: FI1<->FI2 !! C SUCH CHANGE MEANS THAT ANGLES FOUND SHOULD BE REVERSED TOO: C ALPHA1<->ALPHA2 DIMENSION FI1(3),FI2(3) DOUBLE PRECISION R1(3,3),R2(3,3),R3(3,3) DOUBLE PRECISION PSI,THETA,PHI,DEPS CHARACTER*1 MODIS COMMON /POLS/ MODIS DOUBLE PRECISION QUADPI,RAD_TO_DGR PARAMETER (QUADPI = 3.141592653589793238462643383279502884197) PARAMETER (RAD_TO_DGR = (180.0/QUADPI)) PARAMETER (DEPS = 1.0D-7) CALL BLDR(R1,-FI1(3),-FI1(2),-FI1(1)) CALL BLDR(R2,FI2(1),FI2(2),FI2(3)) DO 1 I=1,3 DO 1 J=1,3 R3(I,J)=0.0 DO 1 K=1,3 1 R3(I,J)=R3(I,J)+R2(I,K)*R1(K,J) C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C LIMIT PRECISION DO 5 J=1,3 DO 5 I=1,3 IF(DABS(R3(I,J)).LT.DEPS) R3(I,J)=0.0D0 IF(R3(I,J)-1.0D0.GT.-DEPS) R3(I,J)=1.0D0 IF(R3(I,J)+1.0D0.LT.DEPS) R3(I,J)=-1.0D0 5 CONTINUE IF(R3(3,3).EQ.1.0) THEN THETA=0.0 PSI=0.0 IF(R3(1,1).EQ.0.0) THEN PHI=RAD_TO_DGR*DASIN(R3(1,2)) ELSE PHI=RAD_TO_DGR*DATAN2(R3(1,2),R3(1,1)) ENDIF ELSEIF(R3(3,3).EQ.-1.0) THEN THETA=180.0 PSI=0.0 IF(R3(1,1).EQ.0.0) THEN PHI=RAD_TO_DGR*DASIN(-R3(1,2)) ELSE PHI=RAD_TO_DGR*DATAN2(-R3(1,2),-R3(1,1)) ENDIF ELSE THETA=RAD_TO_DGR*DACOS(R3(3,3)) ST=DSIGN(1.0D0,THETA) IF(R3(3,1).EQ.0.0) THEN IF(ST.NE.DSIGN(1.0D0,R3(3,2))) THEN PHI=270.0 ELSE PHI=90.0 ENDIF ELSE PHI=RAD_TO_DGR*DATAN2(R3(3,2)*ST,R3(3,1)*ST) ENDIF IF(R3(1,3).EQ.0.0) THEN IF(ST.NE.DSIGN(1.0D0,R3(2,3))) THEN PSI=270.0 ELSE PSI=90.0 ENDIF ELSE PSI=RAD_TO_DGR*DATAN2(R3(2,3)*ST,-R3(1,3)*ST) ENDIF ENDIF IF (PSI.LT.0.0) PSI=PSI+360.0 IF (THETA.LT.0.0) THETA=THETA+360.0 IF (PHI.LT.0.0) PHI=PHI+360.0 C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< C IF(DABS(R3(3,1)).LT.DEPS .AND. DABS(R3(3,2)).LT.DEPS) THEN C PRINT *,' ALPHA 1 CANNOT BE CALCULATED' C DID=2.0 C RETURN C ENDIF C IF(DABS(R3(1,3)).LT.DEPS .AND. DABS(R3(2,3)).LT.DEPS) THEN C PRINT *,' ALPHA 2 CANNOT BE CALCULATED' C DID=2.0 C RETURN C ENDIF ALPHA1=90.0+PHI ALPHA2=90.0-PSI IF(ALPHA1.GE.0.0) THEN L1=MOD(INT(NANG*ALPHA1/180.0),2*NANG)+1 ELSE L1=MOD(INT(NANG*(360.0+ALPHA1)/180.0),2*NANG)+1 ENDIF IF(ALPHA2.GE.0.0) THEN L2=MOD(INT(NANG*ALPHA2/180.0),2*NANG)+1 ELSE L2=MOD(INT(NANG*(360.0+ALPHA2)/180.0),2*NANG)+1 ENDIF END FUNCTION EDL1(X1,X2,N) DIMENSION X1(N),X2(N) TDL1=0.0 DO 1 I=1,N 1 TDL1=TDL1+(X1(I)-X2(I))*(X1(I)-X2(I)) EDL1=TDL1 END DOUBLE PRECISION FUNCTION EDL1_(X1,X2,N) DIMENSION X1(N),X2(N) DOUBLE PRECISION TDL1 TDL1=0.0 DO 1 I=1,N 1 TDL1=TDL1+(X1(I)-X2(I))*DBLE(X1(I)-X2(I)) EDL1_=TDL1 END C DOUBLE PRECISION FUNCTION EDL2(X1,X2,N) DIMENSION X1(N),X2(N) DOUBLE PRECISION TDL2 TDL2=0.0 DO 1 I=1,N,2 TDL2=TDL2+(X1(I)-X2(I))*DBLE(X1(I)-X2(I)) 1 TDL2=TDL2+(X1(I+1)-X2(I+1))*DBLE(X1(I+1)+X2(I+1)) EDL2=TDL2 END C DOUBLE PRECISION FUNCTION EDL3(X1,X2,N) DIMENSION X1(N),X2(N) DOUBLE PRECISION TDL3 TDL3=0.0 DO 1 I=1,N,2 TDL3=TDL3+(X1(I)-X2(I))*DBLE(X1(I)-X2(I)) 1 TDL3=TDL3+(X1(I+1)-X2(I+1))*DBLE(X1(I+1)+X2(I+1)) EDL3=TDL3 END C DOUBLE PRECISION FUNCTION ECL1(X1,X2,N) DIMENSION X1(N),X2(N) DOUBLE PRECISION TCL1 TCL1=1.0D0 DO 1 I=1,N 1 TCL1=TCL1-X1(I)*X2(I) ECL1=TCL1 END C DOUBLE PRECISION FUNCTION ECL2(X1,X2,N) DIMENSION X1(N),X2(N) DOUBLE PRECISION TCL2 TCL2=1.0D0 DO 1 I=1,N 1 TCL2=TCL2-X1(N-I+1)*X2(I) ECL2=TCL2 END C DOUBLE PRECISION FUNCTION ECL3(X1,X2,N) DIMENSION X1(N),X2(N) DOUBLE PRECISION TEMP TEMP=1.0D0 DO 1 I=1,N 1 TEMP=TEMP-X1(I)*X2(N-I+1) ECL3=TEMP END C SUBROUTINE NRL(X,M) DIMENSION X(M) DOUBLE PRECISION A,S A=0.0 S=0.0 A=SUM(X) S=DOT_PRODUCT(X,X) A=A/M S=DSQRT(S-M*A*A) X=(X-A)/S END SUBROUTINE SORT2(RA,RB,N) INTEGER RA(N),RB(N),RRA,RRB L=N/2+1 IR=N 10 CONTINUE IF (L.GT.1)THEN L=L-1 RRA=RA(L) RRB=RB(L) ELSE RRA=RA(IR) RRB=RB(IR) RA(IR)=RA(1) RB(IR)=RB(1) IR=IR-1 IF(IR.EQ.1)THEN RA(1)=RRA RB(1)=RRB RETURN ENDIF ENDIF I=L J=L+L 20 IF(J.LE.IR)THEN IF(J.LT.IR)THEN IF(RA(J).LT.RA(J+1))J=J+1 ENDIF IF(RRA.LT.RA(J))THEN RA(I)=RA(J) RB(I)=RB(J) I=J J=J+J ELSE J=IR+1 ENDIF GO TO 20 ENDIF RA(I)=RRA RB(I)=RRB GO TO 10 END SUBROUTINE PREPANG(DM,NANG) DIMENSION DM(2,NANG) DOUBLE PRECISION PHI PHI=4.0D0*DATAN(1.0D0)/NANG DO I=1,NANG DM(1,I)=DCOS(PHI*(I-1)) DM(2,I)=DSIN(PHI*(I-1)) ENDDO C USE NEXT TWO LINES TO GET THE SAME SINOGRAM AS IN _RM 2D_ C DM(1,I)=-DSIN(PHI*(I-1)) C1 DM(2,I)=DCOS(PHI*(I-1)) END SUBROUTINE RDL_P(N,NANG,LENPROJ,RI,INUMBRT,NIMA,PRJ,LENB,LENF) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' DIMENSION PRJ(LENF,2*NANG,NIMA) DIMENSION INUMBRT(NIMA) C AUTOMATIC ARRAYS, PROJECTION IS ASSUMED TO HAVE ODD LENGTH DIMENSION B1(LENPROJ+1),B2(LENPROJ+15),BT(LENPROJ+1) REAL, ALLOCATABLE, DIMENSION(:,:) :: X,PROJ,DM INTEGER, ALLOCATABLE, DIMENSION(:,:) :: IPCUBE LOGICAL BREAK COMMON /F_SPEC/ FINPAT,FINPIC,FINFO,FINDOC,NLET CHARACTER*80 FINPAT,FINPIC,FINFO,FINDOC COMMON /PI/ PI DATA INPIC/99/ DATA NDOUT/89/ PI=4.0*DATAN(1.0D0) C CALCULATE DIMENSIONS FOR NORMAS NSB=-N/2 NSE=-NSB-1+MOD(N,2) NRB=-N/2 NRE=-NRB-1+MOD(N,2) C RADIUS FOR NORMAS IR1=0 IR2=LENPROJ/2 ALLOCATE(DM(2,NANG),STAT=IRTFLG) IF (IRTFLG.NE.0) CALL ERRT(46,'OP, DM',IER) ALLOCATE(X(N,N),STAT=IRTFLG) IF (IRTFLG.NE.0) CALL ERRT(46,'OP, X',IER) ALLOCATE(PROJ(LENPROJ,NANG),STAT=IRTFLG) IF (IRTFLG.NE.0) CALL ERRT(46,'OP, PROJ',IER) CALL PREPANG(DM,NANG) C FIRST CALL JUST FIGURES NN, INUMBRT HERE IS A DUMMY ARGUMENT CALL PSQU_P(N,NN,INUMBRT,RI,.FALSE.) ALLOCATE(IPCUBE(3,NN),STAT=IRTFLG) IF (IRTFLG.NE.0) CALL ERRT(46,'OP, IPCUBE',IER) CALL PSQU_P(N,NN,IPCUBE,RI,.TRUE.) #ifdef SP_LIBFFT LDA=1 CALL SCFFT1DUI(LENPROJ,B2) #endif DO 8 K=1,NIMA C READ ONE IMAGE CALL FILGET(FINPAT,FINPIC,NLET,INUMBRT(K),INTFLAG) MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINPIC,INPIC,'O',IFORM, & LSAM,LROW,NSL,MAXIM,' ',.FALSE.,NF) IF (NF.NE.0) RETURN DO J=1,N CALL REDLIN(INPIC,X(1,J),N,J) ENDDO CLOSE(INPIC) C NORMALIZE UNDER THE MASK CALL NORMAS(X,NSB,NSE,NRB,NRE,IR1,IR2) C CALCULATE LINE PROJECTIONS CALL PRJC2_P(X,N,DM,NANG,PROJ,LENPROJ,IPCUBE,NN) C-------------------------------------------------------------- BREAK=.FALSE. Cc$omp parallel do private(j,b1,bt),shared(break) DO 1 J=1,NANG B1(1:LENPROJ)=PROJ(1:LENPROJ,J) C NORMALIZE THE LINE CALL NRL(B1,LENPROJ) C MIRROR THE LINES DO L=1,LENPROJ BT(L)=B1(LENPROJ-L+1) ENDDO INV=+1 #ifdef SP_LIBFFT CALL SCFFT1DU(INV,LENPROJ,B1,LDA,B2) #else CALL FMR_1(B1,LENPROJ,B2,INV) #endif IF (INV.NE.1) THEN BREAK=.TRUE. ELSE DO L=LENB,LENB+LENF-1 PRJ(L-LENB+1,J,K)=B1(L) ENDDO ENDIF INV=+1 #ifdef SP_LIBFFT CALL SCFFT1DU(INV,LENPROJ,BT,LDA,B2) #else CALL FMR_1(BT,LENPROJ,B2,INV) #endif DO L=LENB,LENB+LENF-1 PRJ(L-LENB+1,J+NANG,K)=BT(L) ENDDO 1 CONTINUE IF (BREAK) THEN CALL ERRT(38,'OP',NE) RETURN ENDIF 8 CONTINUE DEALLOCATE(IPCUBE) DEALLOCATE(PROJ) DEALLOCATE(X) DEALLOCATE(DM) END SUBROUTINE PRJC2_P (SQUARE,N,DM,NANG,PROJ,LENPROJ,IPCUBE,NN) DIMENSION DM(2,NANG) DIMENSION SQUARE(N,N),PROJ(LENPROJ,NANG) INTEGER IPCUBE(3,NN) COMMON /PAR/ LDP,NM,LDPNM PROJ=0.0 CC$OMP PARALLEL DO PRIVATE(I) DO I=1,NANG CALL PRJC11_P(SQUARE,N,DM(1,I),PROJ(1,I),LENPROJ,IPCUBE,NN) ENDDO END SUBROUTINE PRJC11_P(SQUARE,N,DM,PROJ,LENPROJ,IPCUBE,NN) DIMENSION DM(2) DIMENSION SQUARE(N,N),PROJ(LENPROJ) INTEGER IPCUBE(3,NN) COMMON /PAR/ LDP,NM,LDPNM DO 1 I=1,NN IY=IPCUBE(2,I) XB=(IPCUBE(1,I)-LDP)*DM(1)+(IY-LDP)*DM(2) DO 1 J=IPCUBE(1,I),IPCUBE(3,I) IQX=IFIX(XB+FLOAT(LDPNM)) IF(IQX.GT.0 .AND. IQX.LT.LENPROJ) THEN DIPX=XB+LDPNM-IQX PROJ(IQX)=PROJ(IQX)+(1.0-DIPX)*SQUARE(J,IY) PROJ(IQX+1)=PROJ(IQX+1)+DIPX*SQUARE(J,IY) ENDIF 1 XB=XB+DM(1) END SUBROUTINE PSQU_P(N,NN,IPCUBE,RI,FILL) INTEGER IPCUBE(3,*) LOGICAL FIRST,FILL COMMON /PAR/ LDP C IPCUBE: C 1 - IX C 2 - IY C 3 - END OF IX R=RI*RI NN=0 DO 25 I1=1,N T=I1-LDP XX=T*T FIRST=.TRUE. DO 20 I2=1,N T=I2-LDP RC=T*T+XX IF (FIRST) THEN IF (RC.GT.R) GOTO 14 FIRST = .FALSE. NN = NN+1 IF (FILL) THEN IPCUBE(1,NN)=I2 IPCUBE(2,NN)=I1 IPCUBE(3,NN)=I2 ENDIF ELSE IF (FILL) IPCUBE(3,NN)=I2 IF (RC.GT.R) GOTO 16 ENDIF 14 CONTINUE 20 CONTINUE 16 CONTINUE 25 CONTINUE END SUBROUTINE GENAT(ANGTAB,DPSI,NANGT,FILL) DIMENSION ANGTAB(2,*) LOGICAL FILL PARAMETER (QUADPI = 3.141592653589793238462643383279502884197) PARAMETER (DGR_TO_RAD = (QUADPI/180)) C ANGTAB(1 - THETA C ANGTAB(2 - PHI NANGT=0 P1=0.0 P2=359.9 DO THETA=0.,179.999,DPSI IF (THETA.EQ.0.0.OR.THETA.EQ.180.0) THEN DETPHI=360.0 LT=1 ELSE DETPHI=DPSI/SIN(DGR_TO_RAD*THETA) LT=MAX0(INT((P2-P1)/DETPHI)-1,1) DETPHI=(P2-P1)/LT ENDIF C DO 1 PHI=P1,P2,DETPHI DO I=1,LT PHI = P1+(I-1)*DETPHI NANGT = NANGT+1 IF (FILL) THEN ANGTAB(1,NANGT)=THETA ANGTAB(2,NANGT)=PHI ENDIF ENDDO ENDDO END #ifdef SP_IBMSP3 REAL FUNCTION ETIME(TARRAY) DIMENSION TARRAY(2) C CLOCK RESETS AT ICMAX!! CALL SYSTEM_CLOCK(ICOUNT,ICPSEC,ICMAX) ETIME = ICOUNT / ICPSEC RETURN END #endif