/.JANA4S LOGON /SYSFILE SYSLST=C.JANA4S /ER * /EXEC $FOR1 COMOPT MSGLEVEL=E COMOPT PAD COMOPT OPT=2 COMOPT LIST=(SOURCE,XREF) PROGRAM JANA4S COMMON/KLICE/ NA,NS,NF,CSYM,VYH,NC,NOWR,IWQ,YOMIN,YOMAX,BLKOEF, 1 CENTR,NCYKL,SNLMN,SNLMX,KIM,DIF,KIW,LITE,IFSQ COMMON/UNITS/ W,M40,M50 COMMON/REFL/ YO,YC,A,B,DY,WDYQ,NULOVA,SINTHL,RHO,WT,DER(9708), 1 DC(9708),DYP,IQ COMMON/IND/ H(24),K(24),L(24),T(24),HQ(24),KQ(24),LQ(24),HK2(24) 1 ,HL2(24),KL2(24),FSMM(24),EPS(24),EPSM(24),MM,PHFP COMMON/ATOMY/ FF(32,15),FX(15),FY(15),AW(15) COMMON/CM4/ ATOM(100),AI(100),ISF(100),ITF(100),X(3,100), 1 BETA(6,100),SC(6),KMODX(100),UX(3,4,100),UY(3,4,100) 2 ,PHF(100),QCNTM(100),KMODS(100),A0(100),AX(4,100), 3 AY(4,100),KMODB(100),BX(6,4,100),BY(6,4,100), 4 KFX(100),KFB(100) COMMON/SYMM/ IS(2,4,24),TS(3,24),TAU(24),FSM(24) COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10), 1 BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ, 2 EC(12) COMMON/VAZBY/ NVX,NX(10),NVB,NB(10),RMX(3,3,10),RMB(6,6,10), 1 NVAI,NAI(2,10),NAIXB(10),SUMAI(10) COMMON/STALY/ PI2 COMMON/MODUL/ QUI(3),QUIR(3) COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10), 1 UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10), 2 PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10), 3 SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100), 4 NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10), 5 UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10) COMMON/TEMPAR/SQRP(3),RHOM COMMON NPS(9),NAM(9),NPSS(10),NAMM,NP,NZ,TLUM,Q,KI(9708),PS(9708), 1 AM(75000) EQUIVALENCE (NPSS(10),NPSM) LOGICAL CSYM,NULOVA INTEGER W REAL K,L,KQ,LQ,KL2,KP,LP CHARACTER*1 ISPACE,ISTAR,INR,NSPEC,KSPEC CHARACTER*8 MOLNAME,ATOM CHARACTER*2 SYMBOL DATA ISPACE,ISTAR,INR/' ','*','#'/ DIMENSION HH(3),HHP(3),RNUM(6),RZNUM(6),RDEN(6),RZDEN(6),WRNUM(6), 1 WRZNUM(6),WRDEN(6),WRZDEN(6) CALL PRELIM EFPIP(7)=1.5*RR DO 9000IU=1,NCYKL+1 CALL NEWPG IF(IU.EQ.NCYKL+1) GO TO 400 WRITE(W,107) IU,TLUM GO TO 450 400 WRITE(W,110) 450 CALL CONTR(1,0) CALL CONTRM(1,0) IF(NMOL.GT.0) CALL SETAMP DO 500I=1,NPSM 500 PS(I)=0. DO 700I=1,NAMM 700 AM(I)=0. DO 800I=1,6 RNUM(I)=0. RZNUM(I)=0. RDEN(I)=0. RZDEN(I)=0. WRNUM(I)=0. WRZNUM(I)=0. WRDEN(I)=0. 800 WRZDEN(I)=0. NO=0 NZ=0 NVYH=0 1000 READ(91,333) ISENT,IH1,IK1,IL1,IM1,YO,SIGYO,IQ,SINTHL,NXX 333 FORMAT(I1,4I4,2F9.4,I4,F9.6,I4) IF(ISENT.NE.0) GO TO 8000 NO=NO+1 MM=IM1 MMABS=IABS(MM) PHFP=-MMABS*(MMABS-1) IF(KIM.LE.-2) GO TO 1100 IF(KIM.EQ.-1.AND.MM.NE.0) GO TO 1100 IF(MMABS.EQ.KIM) GO TO 1100 GO TO 1000 1100 IF(NXX.GT.NC) GO TO 1000 IF(SINTHL.LT.SNLMN.OR.SINTHL.GT.SNLMX) GO TO 1000 MMABS=MMABS+2 NZ=NZ+1 HH(1)=IH1 HH(2)=IK1 HH(3)=IL1 FMM=MM DO 1110I=1,3 1110 HH(I)=HH(I)+FMM*(QUI(I)+QUIR(I)) POM=1./SQRT(HH(1)**2+HH(2)**2+HH(3)**2) DO 2111I=1,3 EFPIP(I)=POM*HH(I) 2111 EFPIP(I+3)=EFPIP(I) SIGYO=SQRT(SIGYO**2+BLKOEF*YO**2) RHO=SINTHL**2*RHOM NULOVA=NXX.GT.0 NSPEC=ISPACE IF(NULOVA) NSPEC=ISTAR C DO 2500I=1,NS DO 2200J=1,3 2200 HHP(J)=0. DO 2400J=1,3 DO 2300M=1,2 ISP=IS(M,J,I) ISPA=IABS(ISP) IF(ISP) 2230,2300,2250 2230 HHP(ISPA)=HHP(ISPA)-HH(J) GO TO 2300 2250 HHP(ISPA)=HHP(ISPA)+HH(J) 2300 CONTINUE 2400 CONTINUE C EPS(I)=1. IF(IS(1,4,I).LT.0) EPS(I)=-1. IF(IS(1,4,I).EQ.0.AND.IS(2,4,I).LT.0) EPS(I)=-1. FSMM(I)=FSM(I)*FMM EPSM(I)=EPS(I)*FMM HP=HHP(1) KP=HHP(2) LP=HHP(3) H(I)=HP*PI2 K(I)=KP*PI2 L(I)=LP*PI2 T(I)=(TS(1,I)*HH(1)+TS(2,I)*HH(2)+TS(3,I)*HH(3))*PI2 HP=HP*SQRP(1) KP=KP*SQRP(2) LP=LP*SQRP(3) HQ(I)=HP**2 KQ(I)=KP**2 LQ(I)=LP**2 HK2(I)=2.*HP*KP HL2(I)=2.*HP*LP 2500 KL2(I)=2.*KP*LP PT=SINTHL/.05+1. IPT=PT FRACPT=PT-FLOAT(IPT) DO 2800I=1,NF FX(I)=FF(IPT,I)+(FF(IPT+1,I)-FF(IPT,I))*FRACPT+FF(29,I) 2800 FY(I)=FF(31,I)/FX(I) CALL CALC IF(IU.GT.NCYKL) THEN WRITE(80) ISENT,IH1,IK1,IL1,IM1,YO,YC,A,B,SINTHL,NXX GO TO 2990 ENDIF IF(NVAI.EQ.0) GO TO 2819 DO 2818I=1,NVAI J1=NAI(1,I)*92-73 J2=NAI(2,I)*92-73 IF(NAIXB(I).LE.0) GO TO 2810 DER(J1)=DER(J1)-DER(J2) DER(J2)=0. GO TO 2815 2810 DO 2811J=10,18 DER(J1+J)=DER(J1+J)-DER(J2+J) 2811 DER(J2+J)=0. 2815 IF(IABS(NAIXB(I)).EQ.3) GO TO 2818 JK=9 IF(IABS(NAIXB(I)).EQ.2) JK=3 DO 2817J=1,JK DER(J1+J)=DER(J1+J)+DER(J2+J) 2817 DER(J2+J)=0. IF(JK.NE.9) GO TO 2818 DO 12817J=19,91 DER(J1+J)=DER(J1+J)+DER(J2+J) 12817 DER(J2+J)=0. 2818 CONTINUE 2819 IF(NVX.LE.0.AND.NVB.LE.0) GO TO 2900 IF(NVX.LE.0) GO TO 2850 DO 2820I=1,NVX M=NX(I)*92-72 2820 CALL ZMDER(DER(M),RMX(1,1,I),3) 2850 IF(NVB.LE.0) GO TO 2900 DO 2860I=1,NVB M=NB(I)*92-69 2860 CALL ZMDER(DER(M),RMB(1,1,I),6) 2900 IF(NMOL.EQ.0) GO TO 2990 DO 2910I=NPMP,NPMK 2910 DER(I)=0. J1=-73 DO 2960I=1,NA J1=J1+92 IF(KMODX(I).GE.0) GO TO 2960 N=-KMODX(I) M=KMODXM(N) IF(M.LE.0) GO TO 2960 J2=NPMP+49*(N-1) JX=J1+1 JM=J1+19 DO 2920JJ=1,M CALL CULTM(DER(JM),DURDX(1,JJ,N),DER(JX),1,6,3) 2920 JM=JM+6 JM=J2-1 JX=J1+18 JJ=6*M DO 2930J=1,JJ 2930 DER(JM+J)=DER(JM+J)+DER(JX+J) JJ=J1+19 JD=J2-J1+5 DO 2950JX=JJ,JJ+6*(M-1),6 ISP=JX+JD CALL CULTM(DER(JX),DURDR(1,I),DER(ISP),1,3,3) 2950 CALL CULTM(DER(JX+3),DURDR(1,I),DER(ISP+3),1,3,3) JM=J2+48 DER(JM)=DER(JM)+DER(J1+91) 2960 CONTINUE 2990 CALL DSETA(DER(19),NA) IF(NMOL.GT.0) CALL DSETM(DER(NPMP),NMOL) C GO TO (3000,3100,3200),IWQ 3000 WT=1./SIGYO**2 GO TO 3300 3100 WT=1. SIGYO=1. GO TO 3300 3200 WT=1./(YOMIN+YO+YO**2*YOMAX) C 3300 DY=YO-YC DYP=DY YOP=YO**2 IF(IFSQ.NE.2) GO TO 3350 DYP=YOP-YC**2 WT=WT/(4.*YOP) YOP=4.*YOP**2 IF(IU.GT.NCYKL) GO TO 3350 DO 3310I=1,NPMK 3310 DER(I)=DER(I)*2.*YC 3350 WDY=SQRT(WT)*DYP WDYQ=WDY**2 WYOQ=WT*YOP ADY=ABS(DY) AYO=ABS(YO) KSPEC=ISPACE IF(ABS(WDY).LE.VYH) GO TO 4500 KSPEC=INR NVYH=NVYH+1 C 4500 WRZNUM(1)=WRZNUM(1)+WDYQ RZNUM(1)=RZNUM(1)+ADY RZDEN(1)=RZDEN(1)+AYO WRZDEN(1)=WRZDEN(1)+WYOQ IF(NULOVA) GO TO 4700 RNUM(1)=RNUM(1)+ADY RDEN(1)=RDEN(1)+AYO WRNUM(1)=WRNUM(1)+WDYQ WRDEN(1)=WRDEN(1)+WYOQ 4700 WRZNUM(MMABS)=WRZNUM(MMABS)+WDYQ RZNUM(MMABS)=RZNUM(MMABS)+ADY RZDEN(MMABS)=RZDEN(MMABS)+AYO WRZDEN(MMABS)=WRZDEN(MMABS)+WYOQ IF(NULOVA) GO TO 4750 RNUM(MMABS)=RNUM(MMABS)+ADY RDEN(MMABS)=RDEN(MMABS)+AYO WRNUM(MMABS)=WRNUM(MMABS)+WDYQ WRDEN(MMABS)=WRDEN(MMABS)+WYOQ C 4750 CALL CONTR(2,MM) CALL CONTRM(2,MMABS-1) IF(IU.LE.NCYKL) CALL SUMA IF(IABS(NOWR).EQ.2.AND.IU.NE.1.AND.IU.NE.NCYKL+1) GO TO 1000 IF(NOWR) 4800,1000,4900 4800 IF(KSPEC.EQ.ISPACE) GO TO 1000 IF(MOD(NVYH,50).NE.1) GO TO 4950 GO TO 4940 4900 IF(MOD(NZ,50).NE.1) GO TO 4950 4940 CALL NEWPG WRITE(W,100) 4950 WRITE(W,101) IH1,IK1,IL1,IM1,YO,YC,A,B,DY,SIGYO,WDY,NO,NSPEC,KSPEC 1 ,SINTHL,IQ,EXTKOR C GO TO 1000 C 8000 Q=WRZNUM(1) IF(IU.GT.NCYKL) WRITE(80) ISENT,IH1,IK1,IL1,IM1,YO,YC,A,B,SINTHL, 1 NXX CALL NEWPG WRITE(W,108) DO 8500I=1,6 IF(RZDEN(I).EQ.0.) GO TO 8500 IF(I.GE.2) WRITE(W,109) I-2,I-2 R=RNUM(I)/RDEN(I) RZ=RZNUM(I)/RZDEN(I) WRNUM(I)=SQRT(WRNUM(I)) WRDEN(I)=SQRT(WRDEN(I)) WRZNUM(I)=SQRT(WRZNUM(I)) WRZDEN(I)=SQRT(WRZDEN(I)) WR=WRNUM(I)/WRDEN(I) WRZ=WRZNUM(I)/WRZDEN(I) WRITE(W,102) WRITE(W,103) RZNUM(I),RZDEN(I),RZ WRITE(W,104) RNUM(I),RDEN(I),R WRITE(W,105) WRZNUM(I),WRZDEN(I),WRZ WRITE(W,106) WRNUM(I),WRDEN(I),WR 8500 CONTINUE CALL CONTR(3,0) CALL CONTRM(3,0) IF(IU.GT.NCYKL) GO TO 9000 C CALL INVER CALL INTER C 9000 REWIND 91 STOP C 100 FORMAT('0',3X,'H',3X,'K',3X,'L',3X,'M',5X,'YO',9X,'YC',10X,'A',10X 1,'B',9X,'DY',8X,'SIGYO',7X,'WDY',4X,'NREF',4X,'SINTHL',4X,'IQ' 2,6X,'EXT'//) 101 FORMAT(1X,4I4,7F11.4,I5,1X,2A1,F9.5,I5,F9.5) 102 FORMAT(1H0,57X,'NUMERATOR DENOMINATOR R') 103 FORMAT(' R FACTOR INCLUDING UNOBSERVED REFLECTIONS',10X,2F15.3, 1 F10.4) 104 FORMAT(' R FACTOR OMITTING UNOBSERVED REFLECTIONS',10X,2F15.3, 1 F10.4) 105 FORMAT(' WEIGHTED R FACTOR INCLUDING UNOBSERVED REFLECTIONS ', 1 2F15.3,F10.4) 115 FORMAT(' WEIGHTED R FACTOR INCLUDING UNOBSERVED REFLECTIONS ', 1 2F15.3,2F10.4) 106 FORMAT(' WEIGHTED R FACTOR OMITTING UNOBSERVED REFLECTIONS ', 1 2F15.3,F10.4) 107 FORMAT(1X,//////////////////////'0',35X,53('*')/36X,'* REFINEMENT 1CYCLE #',I2,' WITH DUMPING FACTOR : ',F5.2,' *'/36X,53('*')) 108 FORMAT('0TABLE OF R FACTORS FOR ALL REFLECTIONS') 109 FORMAT('0TABLE OF R FACTORS FOR REFLECTIONS WITH SATELLITE INDEX : 1 ', I1,' OR -',I1) 110 FORMAT(1X,//////////////////////'0',51X,21('*')/52X,'* FINAL CALCU 1LATION *'/1X,51X,21('*')) END BLOCK DATA COMMON/UNITS/ W,M40,M50 COMMON/LABEL/ L1,L2(3),L5,L7,L8,L9,L10 COMMON/TITLE/ TITLE(20),NAME(20),PAGE COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10), 1 BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ, 2 EC(12) COMMON/STALY/ PI2 INTEGER W,TITLE,PAGE CHARACTER*12 L1,L2,L5,L7,L8,L9,L10 DATA L1/'OCCUPATION '/,L2/'X ', 'Y ', 1 'Z '/,L5/'U ISO '/, 2 L7/'RHO ISO '/,L8/'G ISO '/,L9/'PHASON '/, 3 L10/'A0 '/, 4 W,M40,M50,PAGE/6,40,50,0/,TORAD/0.017453293/ DATA PI2/6.2831853/ END SUBROUTINE CONTR(ICONT,MM) DIMENSION SSQWDS(8),SSQWDF(9),SINMEZ(8,3),FMEZ(8,3), 1 NLSP(8,3),NLSN(8,3),NLS(8),NLFP(9,3),NLFN(9,3),NLF(9), 2 RNUMSP(8,3),RNUMSN(8,3),RNUMS(8),RDENS(8,3),RFAKS(8), 3 RNUMFP(9,3),RNUMFN(9,3),RNUMF(9),RDENF(9,3),RFAKF(9), 4 SPQWDS(8,3),SNQWDS(8,3),SPQWDF(9,3),SNQWDF(9,3) COMMON/UNITS/ W,M40,M50 COMMON/REFL/ YO,YC,A,B,DY,WDYQ,NULOVA,SINTHL,RHO,WT,DER(9708), 1 DC(9708),DYP,IQ DIMENSION SLIM(3),FLIM(3) LOGICAL NULOVA,CTENO INTEGER W CHARACTER*21 IDENT(3) DATA IDENT/'ALL REFLECTIONS ','MAIN REFLECTIONS ', 1 'SATELLITE REFLECTIONS'/ DATA CTENO/.FALSE./ GO TO (5,15,50),ICONT 5 DO 12J=1,3 DO 10I=1,9 SPQWDF(I,J)=0.0 SNQWDF(I,J)=0.0 RNUMFP(I,J)=0.0 RNUMFN(I,J)=0.0 RDENF(I,J)=0.0 NLFP(I,J)=0 NLFN(I,J)=0 IF(I.EQ.9) GO TO 10 SPQWDS(I,J)=0.0 SNQWDS(I,J)=0.0 RNUMSP(I,J)=0.0 RNUMSN(I,J)=0.0 RDENS(I,J)=0.0 NLSP(I,J)=0 NLSN(I,J)=0 10 CONTINUE 12 CONTINUE IF(CTENO) RETURN READ(M50,75)((SINMEZ(I,J),I=1,8),(FMEZ(I,J),I=1,8),J=1,3) 75 FORMAT (8F9.5/8F9.1) DO 13J=1,3 SLIM(J)=SINMEZ(8,J) 13 FLIM(J)=FMEZ(8,J) CTENO=.TRUE. RETURN 15 DO 30I=1,3 IF(MM.EQ.0.AND.I.EQ.3) GO TO 30 IF(MM.NE.0.AND.I.EQ.2) GO TO 30 IF(YO.GT.FLIM(I).OR.SINTHL.GT.SLIM(I)) RETURN LF=9 IF(NULOVA) GO TO 21 LS=0 17 LS=LS+1 IF(SINTHL.GT.SINMEZ(LS,I)) GO TO 17 LF=0 20 LF=LF+1 IF(YO.GT.FMEZ(LF,I)) GO TO 20 21 IF(DY) 22,24,24 22 RNUMFN(LF,I)=RNUMFN(LF,I)+DY NLFN(LF,I)=NLFN(LF,I)+1 SNQWDF(LF,I)=SNQWDF(LF,I)+WDYQ IF(NULOVA) GO TO 25 RNUMSN(LS,I)=RNUMSN(LS,I)+DY NLSN(LS,I)=NLSN(LS,I)+1 SNQWDS(LS,I)=SNQWDS(LS,I)+WDYQ GO TO 25 24 RNUMFP(LF,I)=RNUMFP(LF,I)+DY NLFP(LF,I)=NLFP(LF,I)+1 SPQWDF(LF,I)=SPQWDF(LF,I)+WDYQ IF(NULOVA) GO TO 25 RNUMSP(LS,I)=RNUMSP(LS,I)+DY NLSP(LS,I)=NLSP(LS,I)+1 SPQWDS(LS,I)=SPQWDS(LS,I)+WDYQ 25 RDENF(LF,I)=RDENF(LF,I)+YO IF(NULOVA) RETURN RDENS(LS,I)=RDENS(LS,I)+YO 30 CONTINUE RETURN 50 DO 70J=1,3 SNUMP=0.0 SNUMN=0.0 FNUMP=0.0 FNUMN=0.0 SDEN=0.0 FDEN=0.0 LSP=0 LSN=0 LFP=0 LFN=0 DO 60I=1,9 RNUMF(I)=RNUMFP(I,J)-RNUMFN(I,J) IF(RDENF(I,J).NE.0.0) RFAKF(I)=100.0*RNUMF(I)/RDENF(I,J) NLF(I)=NLFP(I,J)+NLFN(I,J) SSQWDF(I)=SPQWDF(I,J)+SNQWDF(I,J) IF(NLF(I).NE.0) SSQWDF(I)=SSQWDF(I)/FLOAT(NLF(I)) IF(I.EQ.9) GO TO 60 RNUMS(I)=RNUMSP(I,J)-RNUMSN(I,J) IF(RDENS(I,J).NE.0.0) RFAKS(I)=100.0*RNUMS(I)/RDENS(I,J) SNUMP=SNUMP+RNUMSP(I,J) SNUMN=SNUMN+RNUMSN(I,J) FNUMP=FNUMP+RNUMFP(I,J) FNUMN=FNUMN+RNUMFN(I,J) SDEN=SDEN+RDENS(I,J) FDEN=FDEN+RDENF(I,J) LSP=LSP+NLSP(I,J) LSN=LSN+NLSN(I,J) LFP=LFP+NLFP(I,J) LFN=LFN+NLFN(I,J) NLS(I)=NLSP(I,J)+NLSN(I,J) SSQWDS(I)=SPQWDS(I,J)+SNQWDS(I,J) IF(NLS(I).NE.0) SSQWDS(I)=SSQWDS(I)/FLOAT(NLS(I)) 60 CONTINUE LS=LSP+LSN LF=LFP+LFN IF(LF.LE.0) GO TO 70 CALL NEWPG WRITE(W,54) IDENT(J) WRITE(W,55) WRITE(W,56)(SINMEZ(I,J),I=1,8),(NLSP(I,J),I=1,8),(NLSN(I,J),I=1,8) 1 ,NLS,SSQWDS,(SPQWDS(I,J),I=1,8),(SNQWDS(I,J),I=1,8), 2 (RNUMSP(I,J),I=1,8),(RNUMSN(I,J),I=1,8),RNUMS, 3 (RDENS(I,J),I=1,8),RFAKS WRITE(W,57) WRITE(W,59)(FMEZ(I,J),I=1,8),(NLFP(I,J),I=1,9),(NLFN(I,J),I=1,9) 1 ,NLF,SSQWDF,(SPQWDF(I,J),I=1,9),(SNQWDF(I,J),I=1,9), 2 (RNUMFP(I,J),I=1,9),(RNUMFN(I,J),I=1,9),RNUMF, 3 (RDENF(I,J),I=1,9),RFAKF SNUM=SNUMP-SNUMN FNUM=FNUMP-FNUMN RFS=100.0*SNUM/SDEN RFF=100.0*FNUM/FDEN WRITE(W,58) LSP,LFP,LSN,LFN,LS,LF,SNUMP,FNUMP,SNUMN,FNUMN, 1 SNUM,FNUM,SDEN,FDEN,RFS,RFF 70 CONTINUE 54 FORMAT('0STATISTICS AS A FUNCTION OF SIN(TH)/LAMBDA AND STRUCTURE 1FACTORS FOR ',A21) 55 FORMAT(///1X,'SIN(TH)/LAMBDA') 56 FORMAT(1H+,17X,'LIMITS',5X,8F11.5//18X,'NUMBER +',3X,8I11/25X,'-', 13X,8I11/18X,'TOGETHER',3X,8I11//18X,'AVERAGE WDY',8F11.4/18X, 2'SSQWD +',4X,8F11.4/24X,'-',4X,8F11.4///18X,'NUMERATOR +',8F11.1/ 328X,'-',8F11.1//18X,'TOGETHER',3X,8F11.1/18X,'DENOMINATOR',8F11.1 4//18X,'R FACTOR',3X,8F11.2) 57 FORMAT(///' STRUCT. FACTORS') 59 FORMAT(1H+,17X,'LIMITS',5X,8F11.1,6X,'UNOBS'//18X,'NUMBER +',3X, 19I11/25X,'-',3X,9I11/18X,'TOGETHER',3X,9I11//18X,'AVERAGE WDY', 29F11.4/18X,'SSQWD +',4X,9F11.4/24X,'-',4X,9F11.4///18X,'NUMERATOR 3+',9F11.1/28X,'-',9F11.1//18X,'TOGETHER',3X,9F11.1/18X,'DENOMINATO 4R',9F11.1//18X,'R FACTOR',3X,9F11.2) 58 FORMAT(///32X,16HFINAL CHECK ,15X,14HSIN(TH)/LAMBDA,5X,'STRUCT 1URE FACTORS'//47X,'NUMBER +',9X,I9,13X,I9/54X,1H-,9X,I9,13X,I9/47X 2,'TOGETHER',9X,I9,13X,I9//47X,'NUMERATOR +',5X,F11.1,9X,F13.1/57X, 3'-',5X,F11.1,9X,F13.1//47X,'TOGETHER',8X,F11.1,9X,F13.1/ 447X,'DENOMINATOR',5X,F11.1,9X,F13.1///47X,'R-FACTOR',8X,F11.2,9X, 5F13.2) RETURN END SUBROUTINE CONTRM(ICONT,M) DIMENSION SSQWD(5),SNQWD(5),SPQWD(5),NLP(5),NLN(5),NL(5),RNUMP(5), 1 RNUMN(5),RNUM(5),RDEN(5),RFAK(5) COMMON/UNITS/ W,M40,M50 COMMON/REFL/ YO,YC,A,B,DY,WDYQ,NULOVA,SINTHL,RHO,WT,DER(9708), 1 DC(9708),DYP,IQ INTEGER W LOGICAL NULOVA GO TO (1000,2000,3000),ICONT 1000 DO 1500I=1,5 SPQWD(I)=0. SNQWD(I)=0. RNUMP(I)=0. RNUMN(I)=0. RDEN(I)=0. NLP(I)=0 NLN(I)=0. 1500 CONTINUE RETURN 2000 IF(NULOVA) RETURN IF(DY) 2100,2200,2200 2100 RNUMN(M)=RNUMN(M)+DY NLN(M)=NLN(M)+1 SNQWD(M)=SNQWD(M)+WDYQ GO TO 2300 2200 RNUMP(M)=RNUMP(M)+DY NLP(M)=NLP(M)+1 SPQWD(M)=SPQWD(M)+WDYQ 2300 RDEN(M)=RDEN(M)+YO RETURN 3000 DO 3100I=1,5 RNUM(I)=RNUMP(I)-RNUMN(I) IF(RDEN(I).NE.0.) RFAK(I)=100.*RNUM(I)/RDEN(I) NL(I)=NLN(I)+NLP(I) SSQWD(I)=SNQWD(I)+SPQWD(I) IF(NL(I).NE.0) SSQWD(I)=SSQWD(I)/FLOAT(NL(I)) 3100 CONTINUE CALL NEWPG WRITE(W,100)(I,-I,I=1,4) WRITE(W,101) NLP,NLN,NL,SSQWD,SPQWD,SNQWD,RNUMP,RNUMN,RNUM,RDEN, 1 RFAK RETURN 100 FORMAT('0STATISTICS AS A FUNCTION OF SATELLITE INDEX'//'0SATELLITE 1 INDEX :',22X,'0',10X,4(I2,' OR ',I2,6X)//) 101 FORMAT(19X,'NUMBER +',5I16/26X,'-',5I16/19X,'TOGETHER',5I16//19X, 1 'AVERAGE WDY',5F16.4/19X,'SSQWD +',4X,5F16.4/25X,'-',4X, 2 5F16.4///19X,'NUMERATOR +',5F16.1/29X,'-',5F16.1//19X, 3 'TOGETHER ',5F16.1/19X,'DENOMINATOR',5F16.1//19X,'R FACTO 4R' ,3X,5F16.2) END SUBROUTINE PRELIM COMMON/KLICE/ NA,NS,NF,CSYM,VYH,NC,NOWR,IWQ,YOMIN,YOMAX,BLKOEF, 1 CENTR,NCYKL,SNLMN,SNLMX,KIM,DIF,KIW,LITE,IFSQ COMMON/UNITS/ W,M40,M50 COMMON/CM4/ ATOM(100),AI(100),ISF(100),ITF(100),X(3,100), 1 BETA(6,100),SC(6),KMODX(100),UX(3,4,100),UY(3,4,100) 2 ,PHF(100),QCNTM(100),KMODS(100),A0(100),AX(4,100), 3 AY(4,100),KMODB(100),BX(6,4,100),BY(6,4,100), 4 KFX(100),KFB(100) COMMON/CSM4/ SAI(100),SX(3,100),SB(6,100),SUX(3,4,100), 1 SUY(3,4,100),SPHF(100),SA0(100),SAX(4,100), 2 SAY(4,100),SBX(6,4,100),SBY(6,4,100) COMMON/TITLE/ TITLE(20),NAME(20),PAGE COMMON/ATOMY/ FF(32,15),FX(15),FY(15),AW(15) COMMON/SYMM/ IS(2,4,24),TS(3,24),TAU(24),FSM(24) COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10), 1 BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ, 2 EC(12) COMMON/VAZBY/ NVX,NX(10),NVB,NB(10),RMX(3,3,10),RMB(6,6,10), 1 NVAI,NAI(2,10),NAIXB(10),SUMAI(10) COMMON/MODUL/ QUI(3),QUIR(3) COMMON/STALY/ PI2 COMMON/METRICS/ G(9),GI(9),TR(9),VOL COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10), 1 UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10), 2 PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10), 3 SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100), 4 NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10), 5 UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10) COMMON/COVC/ ICOV COMMON/LABEL/ L1,L2(3),L5,L7,L8,L9,L10 COMMON/TEMPAR/SQRP(3),RHOM COMMON NPS(9),NAM(9),NPSS(10),NAMM,NP,NZ,TLUM,Q,KI(9708),PS(9708), 1 AM(75000) DIMENSION RMT(3,10),URCP(6),RCP(6),PRCP(6) LOGICAL CSYM INTEGER W,TITLE,PAGE CHARACTER*8 ATOM,AT,ATO,MOLNAME CHARACTER*1 IC(7),ICENT,IQS(11),IQUIR,SMB(2),PLUS,MINUS,MEZERA, 1 SMBT(5) CHARACTER*12 L1,L2,L5,L7,L8,L9,L10 EQUIVALENCE (A,RCP(1)),(B,RCP(2)),(C,RCP(3)),(CSA,RCP(4)), 1 (CSB,RCP(5)),(CSG,RCP(6)) DATA IC/'P','A','B','C','I','R','F'/,IQS/'P','A','B','C','L','M', 1 'N','U','V','W','R'/PLUS,MINUS,MEZERA,SMBT/'+','-',' ', 2 '1','S','T','Q','H'/,TPISQ,EPISQ/19.7392088,78.95683521/ C READ(M50,100) NAME READ(M50,129) A,B,C,ALFA,BETT,GAMA,FLAME CSA=COS(TORAD*ALFA) CSB=COS(TORAD*BETT) CSG=COS(TORAD*GAMA) SNB=SQRT(1.-CSB**2) SNG=SQRT(1.-CSG**2) POM=(CSA-CSB*CSG)/SNG TR(1)=A TR(2)=0. TR(3)=0. TR(4)=B*CSG TR(5)=B*SNG TR(6)=0. TR(7)=C*CSB TR(8)=C*POM TR(9)=C*SNB*SQRT(1.-(POM/SNB)**2) VOL=A*B*C*SQRT(1.-CSA**2-CSB**2-CSG**2+2.*CSA*CSB*CSG) VOLUME=VOL**2*12.593E-4/FLAME**3 G(1)=A**2 G(2)=A*B*CSG G(3)=A*C*CSB G(4)=G(2) G(5)=B**2 G(6)=B*C*CSA G(7)=G(3) G(8)=G(6) G(9)=C**2 CALL SMI3(G,GI,CENTR) CALL RECIP(RCP,RCP) DO 250I=1,3 250 PRCP(I)=.25*RCP(I)**2 PRCP(4)=.25*RCP(1)*RCP(2)*RCP(6) PRCP(5)=.25*RCP(1)*RCP(3)*RCP(5) PRCP(6)=.25*RCP(2)*RCP(3)*RCP(4) DO 300I=1,3 300 URCP(I)=TPISQ*RCP(I)**2 URCP(4)=TPISQ*RCP(1)*RCP(2) URCP(5)=TPISQ*RCP(1)*RCP(3) URCP(6)=TPISQ*RCP(2)*RCP(3) RHOM=EPISQ DO 301I=1,3 301 SQRP(I)=SQRT(URCP(I)) C READ(M50,113) QUI READ(M50,101) NS,NI,ICENT,IQUIR CSYM=NI.EQ.1 DO 1000NC=1,7 IF(ICENT.EQ.IC(NC)) GO TO 1100 1000 CONTINUE NC=1 1100 GO TO (1200,1210,1210,1210,1210,1220,1230),NC 1200 CENTR=1. GO TO 1300 1210 CENTR=2. GO TO 1300 1220 CENTR=3. GO TO 1300 1230 CENTR=4. 1300 CENTR=CENTR*FLOAT(3-NI) DO 1310I=1,11 IF(IQUIR.EQ.IQS(I)) GO TO 1320 1310 CONTINUE I=1 1320 IF(I.GE.8) GO TO 1330 DO 1321J=1,3 1321 QUIR(J)=0. IF(I.EQ.1) GO TO 1350 I=I-1 P=.5 IF(I.LT.4) GO TO 1322 I=I-3 P=1. 1322 QUIR(I)=P GO TO 1350 1330 IF(I.EQ.11) GO TO 1340 DO 1335J=1,3 1335 QUIR(J)=.5 I=I-7 QUIR(I)=0. GO TO 1350 1340 QUIR(1)=.333333333 QUIR(2)=.333333333 QUIR(3)=0. C 1350 DO 1360I=1,NS READ(M50,102)(TS(J,I),(IS(K,J,I),K=1,2),J=1,3),SMB, 1 (IS(K,4,I),K=1,2) P=1. J=2 IF(SMB(1).EQ.MEZERA.OR.SMB(1).EQ.PLUS) GO TO 1352 IF(SMB(1).EQ.MINUS) GO TO 1351 J=1 GO TO 1352 1351 P=-1. 1352 DO 1353K=1,5 IF(SMB(J).EQ.SMBT(K)) GO TO 1354 1353 CONTINUE K=1 1354 IF(K.EQ.1) P=0. IF(K.EQ.5) K=6 TAU(I)=P/FLOAT(K) 1360 FSM(I)=PI2*(TAU(I)-QUI(1)*TS(1,I)-QUI(2)*TS(2,I)-QUI(3)*TS(3,I)) C READ(M50,100) READ(M50,100) READ(M50,100) READ(M50,100) READ(M50,101) NF READ(M50,103)((FF(J,I),J=1,32),I=1,NF) READ(M50,103)(AW(I),I=1,NF) CALL CONTR(1,0) CALL NAJDI(M50) READ(M50,100) TITLE CALL NEWPG WRITE(W,108) ICENT IF(CSYM) GO TO 1400 WRITE(W,105) GO TO 1450 1400 WRITE(W,106) 1450 WRITE(W,107) WRITE(W,109)((TS(J,I),(IS(K,J,I),K=1,2),J=1,3),TAU(I), 1 (IS(K,4,I),K=1,2),I=1,NS) WRITE(W,156) QUI WRITE(W,157) QUIR,IQUIR CALL NEWPG WRITE(W,110) WRITE(W,111)(I,(FF(J,I),J=1,32),I=1,NF) READ(M50,112) IWQ,NC,NOWR,ICOV,LITE,IFSQ READ(M50,112) IEXT,ITYPEX,IDISTR,IMOSAIC,IAB READ(M50,112) KIM,KIW IF(KIW.NE.1) GO TO 1470 DO 1460I=1,NF 1460 AW(I)=1. 1470 IF(IEXT.GT.1) ITYPEX=IEXT-1 IF(IEXT.LE.0) IAB=0 IDISTR=MIN0(IDISTR,2) IDISTR=MAX0(IDISTR,1) IMOSAIC=MIN0(IMOSAIC,1) IMOSAIC=MAX0(IMOSAIC,0) READ(M50,113) YOMIN,YOMAX,VYH,BLKOEF,RR,UHMON,SNLMN,SNLMX IF(SNLMX.LE.0.) SNLMX=10. IF(IAB.NE.0) READ(M50,130) AA,BA,DADMI CALL NEWPG IF(IFSQ.NE.2) WRITE(W,180) IF(IFSQ.EQ.2) WRITE(W,181) IWQ=IWQ+1 GO TO (1500,1510,1520),IWQ 1500 WRITE(W,114) BLKOEF BLKOEF=(BLKOEF*.01)**2 GO TO 2000 1510 WRITE(W,115) GO TO 2000 1520 WRITE(W,116) YOMIN,YOMAX YOMIN=2.*YOMIN YOMAX=2./YOMAX 2000 WRITE(W,117) VYH WRITE(W,118) NC IF(KIM.EQ.0) WRITE(W,158) IF(KIM.EQ.-1) WRITE(W,159) IF(KIM.GT.0) WRITE(W,160) KIM,-KIM WRITE(W,150) SNLMN,SNLMX IF(NOWR) 12000,22000,32000 12000 WRITE(W,119) GO TO 32000 22000 WRITE(W,173) GO TO 52000 32000 I=ABS(NOWR) 52000 IF(I.EQ.1) WRITE(W,174) IF(I.EQ.2) WRITE(W,175) IF(ICOV.NE.0) WRITE(W,171) READ(M40,112) NA,NMOL IF(NMOL.LT.0) NMOL=0 IF(NMOL.EQ.0) GO TO 2004 IF(KIW-1) 2001,2002,2003 2001 WRITE(W,165) GO TO 2004 2002 WRITE(W,166) WRITE(W,167) GO TO 2004 2003 WRITE(W,166) WRITE(W,168) 2004 READ(M40,131) SC,(KI(M),M=1,6) READ(M40,131) TO READ(M40,131)(EC(M),M=1,6),(KI(M),M=7,12),(EC(M), 1 M=7,12),(KI(M),M=13,18) L=18 DO 2010I=1,NA K=L+1 L=L+10 READ(M40,121) ATOM(I),ISF(I),KFB(I),KFX(I),ITF(I),KMODX(I), 1 KMODS(I),AI(I),(X(J,I),J=1,3),(BETA(J,I),J=1,6), 2 (KI(M),M=K,L) KMODB(I)=IABS(KMODX(I)/10) M=MOD(IABS(KMODX(I)),10) IF(KMODX(I).GE.0) THEN KMODX(I)=M ELSE KMODX(I)=-M ENDIF QCNTM(I)=-PI2*((QUI(1)+QUIR(1))*X(1,I)+(QUI(2)+QUIR(2))*X(2,I)+ 1 (QUI(3)+QUIR(3))*X(3,I)) K=L+1 L=L+1 A0(I)=1. KI(L)=0 IF(KMODS(I).GT.0) READ(M40,151) A0(I),KI(L) DO 2005J=1,4 K=L+1 L=L+2 IF(J.GT.KMODS(I)) THEN AX(J,I)=0. AY(J,I)=0. KI(K)=0 KI(L)=0 ELSE READ(M40,182) AX(J,I),AY(J,I),KI(K),KI(L) ENDIF 2005 CONTINUE DO 2009J=1,4 K=L+1 L=L+6 IF(J.GT.KMODX(I)) THEN DO 2006M=1,3 UX(M,J,I)=0. 2006 UY(M,J,I)=0. DO 2007M=K,L 2007 KI(M)=0 ELSE READ(M40,131)(UX(M,J,I),M=1,3),(UY(M,J,I),M=1,3),(KI(M),M=K,L) ENDIF 2009 CONTINUE DO 2008J=1,4 K=L+1 L=L+12 IF(J.GT.KMODB(I)) THEN DO 12008M=1,6 BX(M,J,I)=0. 12008 BY(M,J,I)=0. DO 22008M=K,L 22008 KI(M)=0 ELSE READ(M40,131)(BX(M,J,I),M=1,6),(KI(M),M=K,K+5), 1 (BY(M,J,I),M=1,6),(KI(M),M=K+6,L) ENDIF 2008 CONTINUE K=L+1 L=L+1 PHF(I)=0. KI(L)=0 IF(KMODX(I)+KMODS(I)+KMODB(I).GT.0) READ(M40,151) PHF(I),KI(L) 2010 CONTINUE IF(NMOL.EQ.0) GO TO 2017 DO 2016I=1,NMOL READ(M40,163) MOLNAME(I),NPG(I),KMODXM(I),(XM(J,I),J=1,3) NPG(I)=MAX0(1,NPG(I)) KPG(1,I)=1 IF(NPG(I).GT.1) READ(M40,112)(KPG(J,I),J=2,NPG(I)) DO 2014J=1,4 K=L+1 L=L+6 IF(J.LE.KMODXM(I)) GO TO 2013 DO 2011M=1,3 UTX(M,J,I)=0. 2011 UTY(M,J,I)=0. DO 2012M=K,L 2012 KI(M)=0 GO TO 2014 2013 READ(M40,131)(UTX(M,J,I),M=1,3),(UTY(M,J,I),M=1,3),(KI(M),M=K,L) 2014 CONTINUE DO 12014J=1,4 K=L+1 L=L+6 IF(J.LE.KMODXM(I)) GO TO 12013 DO 12011M=1,3 URX(M,J,I)=0. 12011 URY(M,J,I)=0. DO 12012M=K,L 12012 KI(M)=0 GO TO 12014 12013 READ(M40,131)(URX(M,J,I),M=1,3),(URY(M,J,I),M=1,3),(KI(M),M=K,L) 12014 CONTINUE K=L+1 L=L+1 PHFM(I)=0. KI(L)=0 IF(KMODXM(I).GT.0) READ(M40,151) PHFM(I),KI(L) CALL CENTER(I) QCNTMM=-PI2*((QUI(1)+QUIR(1))*XM(1,I)+(QUI(2)+QUIR(2))*XM(2,I)+ 1 (QUI(3)+QUIR(3))*XM(3,I)) DO 62015M=1,NA IF(KMODX(M).NE.-I) GO TO 62015 QCNTM(M)=QCNTMM 62015 CONTINUE 2016 CONTINUE 2017 DO 2019I=1,NA READ(M40,120,END=2020,ERR=2020) SAI(I),(SX(J,I),J=1,3), 1 (SB(J,I),J=1,6) M=KMODS(I) IF(M.LE.0) GO TO 22017 READ(M40,182,END=2020,ERR=2020) SA0(I) DO 12017J=1,M 12017 READ(M40,182,END=2020,ERR=2020) SAX(J,I),SAY(J,I) 22017 DO 2018J=1,4 IF(J.GT.KMODX(I)) GO TO 32017 2018 READ(M40,131,END=2020,ERR=2020)(SUX(K,J,I),K=1,3), 1 (SUY(K,J,I),K=1,3) 32017 DO 12018J=1,4 IF(J.GT.KMODB(I)) GO TO 22018 READ(M40,131,END=2020,ERR=2020)(SBX(K,J,I),K=1,6) 12018 READ(M40,131,END=2020,ERR=2020)(SBY(K,J,I),K=1,6) 22018 IF(KMODS(I)+KMODX(I)+KMODB(I).GT.0) READ(M40,131,END=2020, 1 ERR=2020) SPHF(I) 2019 CONTINUE IF(NMOL.EQ.0) GO TO 2050 DO 32019I=1,NMOL READ(M40,120,END=2020,ERR=2020) K=KMODXM(I) IF(K.LE.0) GO TO 32019 DO 12019J=1,K 12019 READ(M40,131,END=2020,ERR=2020)(SUTX(M,J,I),M=1,3), 1 (SUTY(M,J,I),M=1,3) DO 22019J=1,K 22019 READ(M40,131,END=2020,ERR=2020)(SURX(M,J,I),M=1,3), 1 (SURY(M,J,I),M=1,3) READ(M40,131,END=2020,ERR=2020) SPHFM(I) 32019 CONTINUE GO TO 2050 2020 DO 2030I=1,NA SAI(I)=0. DO 2021J=1,3 2021 SX(J,I)=0. DO 2022J=1,6 2022 SB(J,I)=0. DO 2024K=1,3 DO 2023J=1,4 SUX(K,J,I)=0. 2023 SUY(K,J,I)=0. 2024 CONTINUE DO 2026K=1,6 DO 2025J=1,4 SBX(K,J,I)=0. 2025 SBY(K,J,I)=0. 2026 CONTINUE SPHF(I)=0. 2030 CONTINUE IF(NMOL.EQ.0) GO TO 2050 DO 2040I=1,NMOL DO 2032K=1,3 DO 2031J=1,4 SUTX(K,J,I)=0. SUTY(K,J,I)=0. SURX(K,J,I)=0. 2031 SURY(K,J,I)=0. 2032 CONTINUE SPHFM(I)=0. 2040 CONTINUE 2050 IF(TO.EQ.0.) GO TO 2200 C DO 2100I=1,NA IF(ITF(I).EQ.1) BETA(1,I)=BETA(1,I)+TO 2100 CONTINUE 2200 F=0. C DO 2300I=1,NA J=ISF(I) 2300 F=F+FF(1,J)*AI(I) F=F*FLOAT(NS)*CENTR C WRITE(W,122) F CALL NEWPG WRITE(W,123) IF(LITE.LT.1) GO TO 2350 WRITE(W,176) IF(LITE.EQ.1) WRITE(W,177) GO TO 2360 2350 WRITE(W,178) IF(LITE.LE.-1) WRITE(W,179) 2360 WRITE(W,162) SC,(KI(I),I=1,6) J=1 L=18 C DO 2450I=1,NA J=J+1 IF(MOD(J,50).NE.1) GO TO 2400 CALL NEWPG WRITE(W,123) 2400 N=IABS(KMODX(I))+10*KMODB(I) IF(KMODX(I).LT.0) N=-N WRITE(W,125) ATOM(I),ISF(I),KFB(I),KFX(I),ITF(I),N,KMODS(I),AI(I), 1 (X(M,I),M=1,3) J=J+1 IF(MOD(J,50).NE.1) GO TO 2410 CALL NEWPG WRITE(W,123) 2410 K=L+1 L=L+10 WRITE(W,161) (BETA(M,I),M=1,6),(KI(M),M=K,L) IF(LITE.EQ.0.OR.LITE.GE.2) GO TO 2418 IF(LITE.LT.0) GO TO 2415 IF(ITF(I).NE.2) GO TO 2412 DO 2411M=1,6 DO 12410J=1,KMODB(I) BX(M,J,I)=BX(M,J,I)/URCP(M) 12410 BY(M,J,I)=BY(M,J,I)/URCP(M) 2411 BETA(M,I)=BETA(M,I)/URCP(M) GO TO 2418 2412 BETA(1,I)=BETA(1,I)/EPISQ GO TO 2418 2415 IF(ITF(I).NE.2) GO TO 2417 DO 2416M=1,6 DO 12415J=1,KMODB(I) BX(M,J,I)=BX(M,J,I)*URCP(M) 12415 BY(M,J,I)=BY(M,J,I)*URCP(M) 2416 BETA(M,I)=BETA(M,I)*URCP(M) GO TO 2418 2417 BETA(1,I)=BETA(1,I)*EPISQ 2418 K=L+1 L=L+1 IF(KMODS(I).GT.0) WRITE(W,124) A0(I),KI(L) DO 2419J=1,4 K=L+1 L=L+2 IF(J.GT.KMODS(I)) GO TO 2419 WRITE(W,183) AX(J,I),AY(J,I),KI(K),KI(L) 2419 CONTINUE DO 2430N=1,4 K=L+1 L=L+6 IF(N.GT.KMODX(I)) GO TO 2430 J=J+1 IF(MOD(J,50).EQ.1) THEN CALL NEWPG WRITE(W,123) ENDIF WRITE(W,162)(UX(M,N,I),M=1,3),(UY(M,N,I),M=1,3),(KI(M),M=K,L) 2430 CONTINUE DO 2432N=1,4 K=L+1 L=L+12 IF(N.GT.KMODB(I)) GO TO 2432 J=J+1 IF(MOD(J,50).EQ.1) THEN CALL NEWPG WRITE(W,123) ENDIF WRITE(W,162)(BX(M,N,I),M=1,6),(KI(M),M=K,K+5) J=J+1 IF(MOD(J,50).EQ.1) THEN CALL NEWPG WRITE(W,123) ENDIF WRITE(W,162)(BY(M,N,I),M=1,6),(KI(M),M=K+6,L) 2432 CONTINUE K=L+1 L=L+1 IF(KMODS(I)+KMODX(I)+KMODB(I).LE.0) GO TO 2450 J=J+1 IF(MOD(J,50).EQ.1) THEN CALL NEWPG WRITE(W,123) ENDIF WRITE(W,124) PHF(I),KI(K) 2450 CONTINUE C IF(NMOL.EQ.0) GO TO 2500 DO 2499I=1,NMOL J=J+1 IF(MOD(J,50).NE.1) GO TO 2451 CALL NEWPG WRITE(W,123) 2451 WRITE(W,164) MOLNAME(I),NPG(I),KMODXM(I),(XM(J,I),J=1,3) DO 2460N=1,4 K=L+1 L=L+6 IF(N.GT.KMODXM(I)) GO TO 2460 J=J+1 IF(MOD(J,50).NE.1) GO TO 2455 CALL NEWPG WRITE(W,123) 2455 WRITE(W,162)(UTX(M,N,I),M=1,3),(UTY(M,N,I),M=1,3),(KI(M),M=K,L) 2460 CONTINUE DO 2470N=1,4 K=L+1 L=L+6 IF(N.GT.KMODXM(I)) GO TO 2470 J=J+1 IF(MOD(J,50).NE.1) GO TO 2465 CALL NEWPG WRITE(W,123) 2465 WRITE(W,162)(URX(M,N,I),M=1,3),(URY(M,N,I),M=1,3),(KI(M),M=K,L) 2470 CONTINUE K=L+1 L=L+1 IF(KMODXM(I).LE.0) GO TO 2499 J=J+1 IF(MOD(J,50).NE.1) GO TO 2480 CALL NEWPG WRITE(W,123) 2480 WRITE(W,124) PHFM(I),KI(K) 2499 CONTINUE C 2500 IF(LITE.LT.0.OR.LITE.GT.2) LITE=2 IF(LITE.EQ.1) LITE=0 IF(LITE.EQ.0) GO TO 22500 L5='B'//L5(2:12) RHOM=1. DO 12500M=1,3 12500 SQRP(M)=1. C 22500 DO 72500I=1,NA IF(ITF(I).NE.3) GO TO 62500 IF(LITE.EQ.2) GO TO 42500 DO 32500J=2,3 32500 BETA(J,I)=BETA(1,I) BETA(4,I)=BETA(1,I)*RCP(6) BETA(5,I)=BETA(1,I)*RCP(5) BETA(6,I)=BETA(1,I)*RCP(4) ITF(I)=2 GO TO 62500 42500 POM=BETA(1,I) DO 52500J=1,6 52500 BETA(J,I)=POM*PRCP(J) ITF(I)=2 62500 IF(ITF(I).NE.1.AND.ITF(I).NE.2) ITF(I)=1 72500 CONTINUE C NVX=0 NVB=0 NVAI=0 2501 READ(M50,144) AT,K,ATO DO 2502M=1,NA IF(ATOM(M).EQ.AT) GO TO 2503 2502 CONTINUE GO TO 2550 2503 DO 2504I=1,NA IF(ATOM(I).EQ.ATO) GO TO 2505 2504 CONTINUE GO TO 2512 2505 NVAI=NVAI+1 IF(M.LT.I) GO TO 2506 WRITE(W,170) AT,ATO NVAI=NVAI-1 GO TO 2501 2506 NAI(1,NVAI)=M NAI(2,NVAI)=I NAIXB(NVAI)=K SUMAI(NVAI)=AI(M)+AI(I) IF(K.LE.0) SUMAI(NVAI)=A0(M)+A0(I) GO TO 2501 2512 IF(K.EQ.2) GO TO 2520 NVX=NVX+1 NX(NVX)=M READ(M50,145)((RMX(I,J,NVX),J=1,3),I) DO 2515I=1,3 POM=0. DO 2514J=1,3 2514 POM=POM+RMX(I,J,NVX)*X(J,M) 2515 RMT(I,NVX)=X(I,M)-POM GO TO 2501 2520 NVB=NVB+1 NB(NVB)=M READ(M50,146)((RMB(I,J,NVB),J=1,6),I=1,6) GO TO 2501 2550 IF(NVX.LE.0.AND.NVB.LE.0.AND.NVAI.LE.0) GO TO 2585 CALL NEWPG WRITE(W,147) IF(NVX.LE.0) GO TO 2570 DO 2560M=1,NVX L=NX(M) 2560 WRITE(W,148) ATOM(L),((RMX(I,J,M),J=1,3),RMT(I,M),I=1,3) 2570 IF(NVB.LE.0) GO TO 2580 DO 2571M=1,NVB L=NB(M) 2571 WRITE(W,149) ATOM(L),((RMB(I,J,M),J=1,6),I=1,6) 2580 IF(NVAI.LE.0) GO TO 2585 DO 2584M=1,NVAI J1=NAI(1,M) J2=NAI(2,M) IF(NAIXB(M).LT.0) GO TO 12583 WRITE(W,152) ATOM(J1),ATOM(J2),SUMAI(M) 12580 GO TO (2581,2582,2583),IABS(NAIXB(M)) 2581 WRITE(W,153) GO TO 2584 2582 WRITE(W,154) GO TO 2584 2583 WRITE(W,155) GO TO 2584 12583 WRITE(W,184) ATOM(J1),ATOM(J2) GO TO 12580 2584 CONTINUE C 2585 IF(NMOL.LE.0) GO TO 2599 J=0 DO 2588I=1,NMOL IF(NPG(I).LE.1) GO TO 2588 IF(J.NE.0) GO TO 2587 CALL NEWPG WRITE(W,169) 2587 WRITE(W,172) MOLNAME(I),(KPG(J,I),J=1,NPG(I)) 2588 CONTINUE C 2599 IF(IEXT.LE.0) GO TO 2900 CALL NEWPG WRITE(W,132) RR,UHMON C2THMQ=COS(2.*UHMON*TORAD)**2 IF(IEXT.EQ.1) WRITE(W,133) IF(IEXT.NE.1) WRITE(W,134) GO TO (2600,2610,2620),ITYPEX 2600 WRITE(W,135) GO TO 2630 2610 WRITE(W,136) GO TO 2630 2620 WRITE(W,137) 2630 IF(ITYPEX.EQ.2) GO TO 2640 IF(IDISTR.EQ.1) WRITE(W,138) IF(IDISTR.EQ.2) WRITE(W,139) 2640 IF(IEXT.EQ.1.OR.IEXT.EQ.3) GO TO 2700 IF(IMOSAIC.EQ.0) WRITE(W,140) IF(IMOSAIC.EQ.1) WRITE(W,141) 2700 IF(IAB.EQ.0) GO TO 2800 WRITE(W,142) AA,BA,DADMI 2800 WRITE(W,143) (EC(I),I=1,6),(KI(I),I=7,12),(EC(I),I=7,12), 1 (KI(I),I=13,18) C 2900 NP=NA*92+18 NPMP=NP+1 NPMK=NP+NMOL*49 J=0 NAMM=0 NPSS(1)=0 DO 5000I=1,9 NPSI=0 DO 3000M=1,NPMK IF(KI(M).EQ.I) NPSI=NPSI+1 3000 CONTINUE IF(NPSI.LE.200) GO TO 4100 WRITE(W,126) I,NPSI STOP 4100 NPS(I)=NPSI J=J+NPSI NPSS(I+1)=J NAM(I)=NPSI*(NPSI+1)/2 NAMM=NAMM+NAM(I) 5000 CONTINUE IF(NAMM.LE.75000) GO TO 7000 WRITE(W,127) STOP C 7000 CALL NEWPG WRITE(W,128) (I,I=1,9),NPS,J,NAM,NAMM READ(M50,104) TLUM,NCYKL IF(TLUM.LE.0.0.OR.TLUM.GT.1.0) TLUM=1.0 READ(M50,113) DIF RETURN 100 FORMAT(20A4) 101 FORMAT(2I5,4X,A1,4X,A1) 102 FORMAT(3(F11.6,2I2),3X,2A1,1X,2I2) 103 FORMAT(8F9.3) 104 FORMAT(F5.2,I5) 105 FORMAT('0NON-CENTROSYMMETRIC SUPER-SPACE GROUP') 106 FORMAT('0CENTROSYMMETRIC SUPER-SPACE GROUP') 107 FORMAT('0SYMMETRY CARDS'//) 108 FORMAT('0CENTERING OF CELL : ',A1) 109 FORMAT(4(F20.6,2I2)) 110 FORMAT('0ATOMIC SCATTERING TABLES :') 111 FORMAT('0',I5,5X,8F9.3/11X,8F9.3 1 /11X,8F9.3 2 /11X,8F9.3) 112 FORMAT(11I5) 113 FORMAT(8F10.5) 114 FORMAT('0WEIGHT 1/SIGMA**2 COEFFICIENT OF UNSTABILITY IS ',F8.2) 115 FORMAT('0UNIT WEIGHT') 116 FORMAT('0WEIGHT AFTER CRUIKSHANCK WITH YOMIN =',F7.2,' YOMAX =', 1 F9.2) 117 FORMAT('0REFLECTIONS WITH ABS(YO-YC)>',F5.2,'*SIGYO WILL BE SYMBOL 1IZED IN OUTPUT BY #') 118 FORMAT('0REFLECTIONS WITH NXX>',I2,' WILL BE OMITTED') 119 FORMAT('0ONLY REFLECTIONS SYMBOLISED BY # (SEE ABOVE) WILL BE PRIN 1TED') 120 FORMAT(18X,4F9.6/6F9.6) 121 FORMAT(A8,I3,3I1,I3,I1,4F9.6/6F9.6,6X,10I1) 122 FORMAT('0F000 =',F8.0) 123 FORMAT('0INPUT PARAMETERS AND THEIR KEYS FOR REFINEMENT'//) 124 FORMAT(1X,F9.6,48X,'=>',I1,'<=') 125 FORMAT(1X,A8,I3,3I1,I3,I1,4F9.6) 126 FORMAT('0ERROR - NUMBER PARAMETERS TO BE REFINED IN ',I2,'.BLOCK I 1S :',I5,' THAT IS GREATER THAN 200') 127 FORMAT('0NUMBER OF ELEMENTS OF MATRIX IS',I10,' THAT IS GREATER TH 1AN 75000') 128 FORMAT('0BLOCK :',9I8,' TOGETHER'/'0NPS :',9I8,I10/' NAM :', 19I8,I10) 129 FORMAT(3F10.3,3F10.2,F10.5) 130 FORMAT(10F7.4) 131 FORMAT(6F9.6,6X,6I1) 132 FORMAT('0INFORMATION ON EXTINCTION CORRECTION :'/39('=')//'0RADIUS 1 OF SHERICAL SAMPLE :',F9.4,' CM'/'0MONOCHROMATOR ANGLE :',F8.2, 2' DEGREES') 133 FORMAT('0ISOTROPIC CORRECTION') 134 FORMAT('0ANISOTROPIC CORRECTION') 135 FORMAT('0EXTINCTION TYP I') 136 FORMAT('0EXTINCTION TYP II') 137 FORMAT('0GENERAL TYP OF EXTINCTION') 138 FORMAT('0GAUSSIAN DISTRIBUTION') 139 FORMAT('0LORENTZIAN DISTRIBUTION') 140 FORMAT('0ANISOTROPY OF MOSAIC SPREAD ACCORDING COPPENS AND HAMILTO 1N') 141 FORMAT('0ANISOTROPY OF MOSAIC SPREAD ACCORDING NELMES AND THORNLEY 1') 142 FORMAT('0MODIFICATION FOR ABSORPTION'/'0A(SIN(TH)) STEP 0.1 IN SIN 1(TH) : ',10F8.2 /' B(SIN(TH)) STEP 0.1 IN SIN 2(TH) : ',10F8.2 /' DADMI(TH) STEP 5 DEGREES 3 IN TH : ',10F8.4/45X,9F8.4) 143 FORMAT('0LIST OF PARAMETERS OF EXTINCTION AND THEIR KEYS OF REFINE 1MENT' /'0RHO : ',6F9.6,5X,'=>',6I1,'<='/ 2 '0G(I,J) : ',6F9.6,5X,'=>',6I1,'<=') 144 FORMAT(A8,I2,2X,A8) 145 FORMAT(3F5.0) 146 FORMAT(6F5.0) 147 FORMAT('0RESTRINCTIONS :'/16('=')) 148 FORMAT('0RESTRINCTION ON COORDINATES OF ATOM :',A8/ 1 '0X = ',F5.1,'*X +',F5.1,'*Y +',F5.1,'*Z +',F10.5/ 2 ' Y = ',F5.1,'*X +',F5.1,'*Y +',F5.1,'*Z +',F10.5/ 3 ' Z = ',F5.1,'*X +',F5.1,'*Y +',F5.1,'*Z +',F10.5) 149 FORMAT('0RESTRINCTION ON TEMPERATURE PARAMETERS OF ATOM :',A8/ 1 '0B(1,1) = ',F5.1,'*B(1,1) +',F5.1,'*B(2,2) +',F5.1,'*B(3,3 2) +',F5.1,'*B(1,2) +',F5.1,'*B(1,3) +',F5.1,'*B(2,3)'/ 1 ' B(2,2) = ',F5.1,'*B(1,1) +',F5.1,'*B(2,2) +',F5.1,'*B(3,3 2) +',F5.1,'*B(1,2) +',F5.1,'*B(1,3) +',F5.1,'*B(2,3)'/ 1 ' B(3,3) = ',F5.1,'*B(1,1) +',F5.1,'*B(2,2) +',F5.1,'*B(3,3 2) +',F5.1,'*B(1,2) +',F5.1,'*B(1,3) +',F5.1,'*B(2,3)'/ 1 ' B(1,2) = ',F5.1,'*B(1,1) +',F5.1,'*B(2,2) +',F5.1,'*B(3,3 2) +',F5.1,'*B(1,2) +',F5.1,'*B(1,3) +',F5.1,'*B(2,3)'/ 1 ' B(1,3) = ',F5.1,'*B(1,1) +',F5.1,'*B(2,2) +',F5.1,'*B(3,3 2) +',F5.1,'*B(1,2) +',F5.1,'*B(1,3) +',F5.1,'*B(2,3)'/ 1 ' B(2,3) = ',F5.1,'*B(1,1) +',F5.1,'*B(2,2) +',F5.1,'*B(3,3 2) +',F5.1,'*B(1,2) +',F5.1,'*B(1,3) +',F5.1,'*B(2,3)') 150 FORMAT('0SIN(TH)/LAMBDA LIMITS FOR ACCEPTANCE OF REFLECTION', 1 2F10.5) 151 FORMAT(F9.6,51X,I1) 152 FORMAT('0OCCUPATIONAL RESTRINCTION OF ATOMS : ',A8,' AND : ',A8, 1 ' SUMMA OF OCCUPATIONS WILL BE : ',F9.6) 153 FORMAT(38X,'COORDINATES ,TEMPERATURE AND POSITIONAL MODULATION PAR 1AMETERS OF THESE ATOMS WILL BE IDENTICAL') 154 FORMAT(38X,'COORDINATE OF THESE ATOMS WILL BE IDENTICAL') 155 FORMAT(38X,'WITHOUT RESTRINCTION ON COORDINATES, TEMPERATURE AND M 1ODULATION PARAMETERS OF THESE ATOMS') 156 FORMAT('0IRRATIONAL PART OF MODULATION VECTOR Q : ',3F9.5) 157 FORMAT('0 RATIONAL PART OF MODULATION VECTOR Q : ',3F9.5,' DERIVE 1D FROM SYMBOL : ',A1) 158 FORMAT('0ONLY MAIN REFLECTIONS WILL BE ACCEPTED INTO CALCULATION') 159 FORMAT('0ONLY SATELLITE REFLECTIONS WILL BE ACCEPTED INTO CALCULAT 1ION') 160 FORMAT('0ONLY REFLECTIONS WITH SATELLITE INDEX',I2,' OR',I2,'WILL 1ACCEPTED INTO CALCULATION') 161 FORMAT(1X,6F9.6,3X,'=>',10I1,'<=') 162 FORMAT(1X,6F9.6,3X,'=>', 6I1,'<=') 163 FORMAT(A8,3X,2I3,10X,3F9.6) 164 FORMAT(1X,A8,3X,2I3,10X,3F9.6) 165 FORMAT('0MOLECULAR CENTERS WILL NOT BE CHANGED') 166 FORMAT('0MOLECULAR CENTERS WILL BE CHANGED TO BE') 167 FORMAT('+',40X,'A GEOMETRICAL CENTER OF ATOMS') 168 FORMAT('+',40X,'A CENTER OF MASS OF ATOMS') 169 FORMAT('0INFORMATION ON A CRYSTALLOGRAPHIC POINT SYMMETRY OF RIGID 1 GROUPS (MOLECULES)'/'0ELEMENTS OF POINT SYMMETRY ARE GIVEN BY NUM 2BER OF SYMMETRY CARD'/' NEGATIVE NUMBER MEANS THE ELEMENT HAS TO B 3E COMBINED WITH INVERSION CENTER'//) 170 FORMAT('OERROR ATOM : ',A8,' IS PLACED BEHIND ATOM : ',A8,' RESTRI 1NTION WILL NOT BE USED') 171 FORMAT('0ELEMENTS OF COVARIATION MATRIX LARGER THAN 0.9 WILL BE PR 1INTED') 172 FORMAT(' GROUP : ',A8,' HAS ELEMENTS :',24I4) 173 FORMAT('0PRINT OF REFLECTIONS SUPRESSED') 174 FORMAT('0PRINT OF REFLECTIONS IN EVERY CYCLE OF REFINEMENT') 175 FORMAT('0PRINT OF REFLECTIONS BEFORE FIRST AND AFTER LAST CYCLE OF 1 REFINEMENT') 176 FORMAT('0INPUT TEMPERATURE PARAMETERS ARE BETA') 177 FORMAT(' AND WILL BE CHANGED TO BE U' 1) 178 FORMAT('0INPUT TEMPERATURE PARAMETERS ARE U') 179 FORMAT(' AND WILL BE CHANGED TO BE BE 1TA') 180 FORMAT('0REFINEMENT BASED ON F') 181 FORMAT('0REFINEMENT BASED ON F**2') 182 FORMAT(2F9.6,42X,2I1) 183 FORMAT(1X,2F9.6,39X,'=>',2I1,'<=') 184 FORMAT('0OCCUPATIONAL MODULATION OF ATOM : ',A8,' AND ATOM : ',A8, 1 ' HAS TO BE OPPOSITE') END SUBROUTINE NAJDI(M) CHARACTER*8 NZ1,NZ2 DATA NZ1/'JANA '/ 1000 READ(M,100) NZ2 100 FORMAT(A8) IF(NZ2.NE.NZ1) GO TO 1000 RETURN END SUBROUTINE NEWPG COMMON/UNITS/ W,M40,M50 COMMON/TITLE/ TITLE(20),NAME(20),PAGE INTEGER TITLE,PAGE,W PAGE=PAGE+1 WRITE(W,100) PAGE,NAME,TITLE 100 FORMAT('1REFINEMENT OF MODULATED STRUCTURE',75X,'PAGE =',I3/ 1 ' STRUCTURE : ',20A4/13X,20A4/) RETURN END SUBROUTINE CALC COMMON/KLICE/ NA,NS,NF,CSYM,VYH,NC,NOWR,IWQ,YOMIN,YOMAX,BLKOEF, 1 CENTR,NCYKL,SNLMN,SNLMX,KIM,DIF,KIW,LITE,IFSQ COMMON/REFL/ YO,YC,A,B,DY,WDYQ,NULOVA,SINTHL,RHO,WT,DER(9708), 1 DC(9708),DYP,IQ COMMON/CM4/ ATOM(100),AI(100),ISF(100),ITF(100),X(3,100), 1 BETA(6,100),SC(6),KMODX(100),UX(3,4,100),UY(3,4,100) 2 ,PHF(100),QCNTM(100),KMODS(100),A0(100),AX(4,100), 3 AY(4,100),KMODB(100),BX(6,4,100),BY(6,4,100), 4 KFX(100),KFB(100) COMMON/IND/ H(24),K(24),L(24),T(24),HQ(24),KQ(24),LQ(24),HK2(24) 1 ,HL2(24),KL2(24),FSMM(24),EPS(24),EPSM(24),MM,PHFP COMMON/ATOMY/ FF(32,15),FX(15),FY(15),AW(15) COMMON/SYMM/ IS(2,4,24),TS(3,24),TAU(24),FSM(24) COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10), 1 BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ, 2 EC(12) COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10), 1 UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10), 2 PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10), 3 SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100), 4 NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10), 5 UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10) COMMON/STALY/ PI2 LOGICAL CSYM,NULOVA,IZO,HARM,EPSL,MAIN DIMENSION DADUXX(4),DADUXY(4),DADUXZ(4),DADUYX(4),DADUYY(4), 1 DADUYZ(4),DBDUXX(4),DBDUXY(4),DBDUXZ(4),DBDUYX(4), 2 DBDUYY(4),DBDUYZ(4),UXX(4),UXY(4),UXZ(4),UYX(4),UYY(4), 3 UYZ(4),BESP(9,4),DBFDU1(9,4),DBFDU2(9,4),DCHIDU1(4), 4 DCHIDU2(4),DSDU1(4),DSDU2(4),CHI(4),DADU1(4),DADU2(4), 5 DBDU1(4),DBDU2(4),DSDCHI(4),SB(5),DSBDAX(4),DSBDAY(4), 6 KSI(4),DKSIDAX(4),DKSIDAY(4),DADAX(4),DADAY(4),DBDAX(4) 7 ,DBDAY(4) DIMENSION DADBX1(4),DADBX2(4),DADBX3(4),DADBX4(4),DADBX5(4), 1 DBDBX1(4),DBDBX2(4),DBDBX3(4),DBDBX4(4),DBDBX5(4), 2 DADBY1(4),DADBY2(4),DADBY3(4),DADBY4(4),DADBY5(4), 3 DBDBY1(4),DBDBY2(4),DBDBY3(4),DBDBY4(4),DBDBY5(4), 4 DADBX6(4),DBDBX6(4),DADBY6(4),DBDBY6(4),BX1(4),BX2(4), 5 BX3(4),BX4(4),BX5(4),BX6(4),BY1(4),BY2(4),BY3(4),BY4(4), 6 BY5(4),BY6(4),BESB(5,4),DBFDB1(5,4),DBFDB2(5,4), 7 DETADB1(4),DETADB2(4),DSDB1(4),DSDB2(4),ETA(4),DADB1(4), 8 DADB2(4),DBDB1(4),DBDB2(4),DSDETA(4),DSDB3(4),DSDU3(4), 9 DBFDB3(5,4),DBFDU3(5,4) CHARACTER*8 ATOM,MOLNAME REAL K,L,KQ,LQ,KL2,KJ,LJ,KQJ,LQJ,KL2J,KSI C A=0. B=0. AANOM=0. BANOM=0. M=19 DO 500I=1,6 500 DER(I)=0. EXTKOR=1. MAIN=MM.EQ.0 DO 6000I=1,NA AII=AI(I) QCNTMI=QCNTM(I) KMODXI=KMODX(I) KFXI=KFX(I) KMODBI=KMODB(I) KFBI=KFB(I) KMODSI=KMODS(I) IF(KMODXI.LT.0) KMODXI=KMODXM(-KMODXI) HARM=KMODXI.LE.0.AND.KMODSI.LE.0.AND.KMODBI.LE.0 IF(AII.NE.0..AND.(.NOT.HARM.OR.MAIN)) GO TO 910 DO 900J=M,M+91 DER(J)=0. 900 DC(J)=0. GO TO 6000 910 XI=X(1,I) YI=X(2,I) ZI=X(3,I) J=ISF(I) FXI=FX(J) FYI=FY(J) IZO=ITF(I).EQ.1 B1=BETA(1,I) EXPIJ=1. IF(IZO) GO TO 1001 B2=BETA(2,I) B3=BETA(3,I) B4=BETA(4,I) B5=BETA(5,I) B6=BETA(6,I) 1001 IF(HARM) GO TO 1003 SB(1)=A0(I) DO 8J=1,4 IF(J.LE.KMODSI) GO TO 6 5 SB(J+1)=0. KSI(J)=0. DSBDAX(J)=0. DSBDAY(J)=0. GO TO 8 6 UQR=AX(J,I)**2+AY(J,I)**2 IF(UQR.LE.0) GO TO 5 U=SQRT(UQR) UQR=1./UQR POM=U*.5 SB(J+1)=POM DSBDAX(J)=AX(J,I)/U*.5 DSBDAY(J)=AY(J,I)/U*.5 KSI(J)=ATAN2(-AX(J,I),AY(J,I)) DKSIDAX(J)=-AY(J,I)*UQR DKSIDAY(J)= AX(J,I)*UQR 8 CONTINUE M1P=0 M1K=0 M2P=0 M2K=0 M3P=0 M3K=0 M4P=0 M4K=0 M5P=0 M5K=0 M6P=0 M6K=0 M7P=0 M7K=0 M8P=0 M8K=0 M9P=-KMODSI M9K= KMODSI IF(KMODXI.EQ.0) GO TO 10 M1P=-4 M1K= 4 IF(KMODXI.EQ.1) GO TO 10 M2P=-4 M2K= 4 IF(KMODXI.EQ.2) GO TO 10 M3P=-4 M3K= 4 IF(KMODXI.EQ.3) GO TO 10 M4P=-4 M4K= 4 10 IF(KMODBI.EQ.0) GO TO 15 M5P=-4 M5K= 4 IF(KMODBI.EQ.1) GO TO 15 M6P=-4 M6K= 4 IF(KMODBI.EQ.2) GO TO 15 M7P=-4 M7K= 4 IF(KMODBI.EQ.3) GO TO 15 M8P=-4 M8K= 4 15 LP=0 LK=4 IF(KMODXI+KMODSI+KMODBI.NE.1) GO TO 20 LP=IABS(MM) LK=IABS(MM) 20 DO 1002J=1,KMODXI UXX(J)=UX(1,J,I) UXY(J)=UX(2,J,I) UXZ(J)=UX(3,J,I) UYX(J)=UY(1,J,I) UYY(J)=UY(2,J,I) 1002 UYZ(J)=UY(3,J,I) DO 11002J=1,KMODBI BX1(J)=BX(1,J,I) BX2(J)=BX(2,J,I) BX3(J)=BX(3,J,I) BX4(J)=BX(4,J,I) BX5(J)=BX(5,J,I) BX6(J)=BX(6,J,I) BY1(J)=BY(1,J,I) BY2(J)=BY(2,J,I) BY3(J)=BY(3,J,I) BY4(J)=BY(4,J,I) BY5(J)=BY(5,J,I) 11002 BY6(J)=BY(6,J,I) PHFI=PHF(I) 1003 DADFII=0. IF(CSYM) GO TO 1005 DBDFII=0. 1005 DADXI=0. DADYI=0. DADZI=0. IF(IZO) GO TO 1010 DADB1I=0. DADB2I=0. DADB3I=0. DADB4I=0. DADB5I=0. DADB6I=0. 1010 IF(CSYM) GO TO 1015 DBDXI=0. DBDYI=0. DBDZI=0. IF(IZO) GO TO 1015 DBDB1I=0. DBDB2I=0. DBDB3I=0. DBDB4I=0. DBDB5I=0. DBDB6I=0. 1015 DADA0=0. DBDA0=0. IF(KMODSI.LE.0) GO TO 1020 DO 1017J=1,KMODSI DADAX(J)=0. DADAY(J)=0. IF(CSYM) GO TO 1017 DBDAX(J)=0. DBDAY(J)=0. 1017 CONTINUE 1020 IF(KMODXI.LE.0) GO TO 1022 DO 1021J=1,KMODXI DADUXX(J)=0. DADUXY(J)=0. DADUXZ(J)=0. DADUYX(J)=0. DADUYY(J)=0. DADUYZ(J)=0. IF(CSYM) GO TO 1021 DBDUXX(J)=0. DBDUXY(J)=0. DBDUXZ(J)=0. DBDUYX(J)=0. DBDUYY(J)=0. 1021 DBDUYZ(J)=0. 1022 DO 1023J=1,KMODBI DADBX1(J)=0. DADBX2(J)=0. DADBX3(J)=0. DADBX4(J)=0. DADBX5(J)=0. DADBX6(J)=0. DADBY1(J)=0. DADBY2(J)=0. DADBY3(J)=0. DADBY4(J)=0. DADBY5(J)=0. DADBY6(J)=0. IF(CSYM) GO TO 1023 DBDBX1(J)=0. DBDBX2(J)=0. DBDBX3(J)=0. DBDBX4(J)=0. DBDBX5(J)=0. DBDBX6(J)=0. DBDBY1(J)=0. DBDBY2(J)=0. DBDBY3(J)=0. DBDBY4(J)=0. DBDBY5(J)=0. 1023 DBDBY6(J)=0. 1025 DADPHF=0. DBDPHF=0. DO 5000J=1,NS HJ =H(J) KJ =K(J) LJ =L(J) EPSJ=EPS(J) EPSL=EPSJ.LE.0. IF(IZO) GO TO 2001 HQJ =HQ(J) KQJ =KQ(J) LQJ =LQ(J) HK2J=HK2(J) HL2J=HL2(J) KL2J=KL2(J) 2001 FI = T(J)+HJ*XI+KJ*YI+LJ*ZI IF(HARM) GO TO 2070 DO 2020KK=1,KMODXI U1=UXX(KK)*HJ+UXY(KK)*KJ+UXZ(KK)*LJ IF(KK.NE.KMODXI.OR.KFXI.EQ.0) THEN U2=UYX(KK)*HJ+UYY(KK)*KJ+UYZ(KK)*LJ UQR=U1**2+U2**2 IF(UQR.LE.0.) GO TO 2006 U=SQRT(UQR) RU=1./U UQR=1./UQR CHI(KK)=ATAN2(U2,U1) DO 2002LL=LP,LK POM=BESJ(U,LL,DIF,0) BESP(LL+5,KK)=POM IF(LL.EQ.0) GO TO 2002 IF(MOD(LL,2).EQ.1) POM=-POM BESP(-LL+5,KK)=POM 2002 CONTINUE POMB=BESJ(U,LK+1,DIF,0)*RU DO 2005LL=LK,LP,-1 POMA=BESP(LL+5,KK)*RU POM=-POMB IF(LL.EQ.0) GO TO 2004 POM=POM+FLOAT(LL)*RU*POMA POMU1=POM*U1 POMU2=POM*U2 DBFDU1(LL+5,KK)=POMU1 DBFDU2(LL+5,KK)=POMU2 IF(MOD(LL,2).EQ.0) GO TO 2003 POMU1=-POMU1 POMU2=-POMU2 2003 DBFDU1(-LL+5,KK)=POMU1 DBFDU2(-LL+5,KK)=POMU2 GO TO 2005 2004 DBFDU1(5,KK)=POM*U1 DBFDU2(5,KK)=POM*U2 2005 POMB=POMA C DCHIDU1(KK)=-UQR*U2 DCHIDU2(KK)= UQR*U1 GO TO 2020 2006 DO 2010LL=LP,LK IF(LL.EQ.0) THEN BESP(LL+5,KK)=1. DBFDU1(LL+5,KK)=0. DBFDU2(LL+5,KK)=0. ELSE IF(LL.EQ.1) THEN BESP( LL+5,KK)=0. BESP(-LL+5,KK)=0. DBFDU1( LL+5,KK)= .5 DBFDU1(-LL+5,KK)=-.5 DBFDU2( LL+5,KK)= .5 DBFDU2(-LL+5,KK)=-.5 ELSE BESP( LL+5,KK)=0. BESP(-LL+5,KK)=0. DBFDU1( LL+5,KK)=0. DBFDU1(-LL+5,KK)=0. DBFDU2( LL+5,KK)=0. DBFDU2(-LL+5,KK)=0. ENDIF ENDIF 2010 CONTINUE CHI(KK)=0. DCHIDU1(KK)=0. DCHIDU2(KK)=0. ELSE CHI(KK)=PI2*UYX(KK) DO 2015LL=LP,LK CALL XMOD(BESP(LL+5,KK),DBFDU1(LL+5,KK),DBFDU2(LL+5,KK), 1 DBFDU3(LL+5,KK),U1,UYY(KK),UYZ(KK),LL,KFXI) CALL XMOD(BESP(-LL+5,KK),DBFDU1(-LL+5,KK),DBFDU2(-LL+5,KK), 1 DBFDU3(-LL+5,KK),U1,UYY(KK),UYZ(KK),-LL,KFXI) 2015 CONTINUE DCHIDU1(KK)=PI2 ENDIF 2020 CONTINUE DO 12020KK=1,KMODBI BB1=-HQJ *BX1(KK)-KQJ *BX2(KK)-LQJ*BX3(KK)-HK2J*BX4(KK) 1 -HL2J*BX5(KK)-KL2J*BX6(KK) IF(KK.NE.KMODBI.OR.KFBI.EQ.0) THEN BB2=-HQJ *BY1(KK)-KQJ *BY2(KK)-LQJ*BY3(KK)-HK2J*BY4(KK) 1 -HL2J*BY5(KK)-KL2J*BY6(KK) UQR=BB1**2+BB2**2 IF(UQR.LE.0.) GO TO 12006 U=SQRT(UQR) RU=1./U UQR=1./UQR ETA(KK)=ATAN2(BB2,BB1) DO 12002LL=LP,LK 12002 BESB(LL+1,KK)=BESJ(U,LL,DIF,1) POMB=BESJ(U,LK+1,DIF,1)*RU DO 12005LL=LK,LP,-1 POMA=BESB(LL+1,KK)*RU POM=POMB IF(LL.NE.0) POM=POM+FLOAT(LL)*RU*POMA DBFDB1(LL+1,KK)=POM*BB1 DBFDB2(LL+1,KK)=POM*BB2 12005 POMB=POMA DETADB1(KK)=-UQR*BB2 DETADB2(KK)= UQR*BB1 GO TO 12020 12006 DO 12010LL=LP,LK IF(LL.EQ.0) THEN BESB(1,KK)=1. DBFDB1(1,KK)=0. DBFDB2(1,KK)=0. ELSE IF(LL.EQ.1) THEN BESB(2,KK)=0. DBFDB1(2,KK)=.5 DBFDB2(2,KK)=.5 ELSE BESB(LL+1,KK)=0. DBFDB1(LL+1,KK)=0. DBFDB2(LL+1,KK)=0. ENDIF ENDIF 12010 CONTINUE ELSE ETA(KK)=PI2*BY1(KK) DO 22002LL=LP,LK CALL BMOD(BESB(LL+1,KK),DBFDB1(LL+1,KK),DBFDB2(LL+1,KK), 1 DBFDB3(LL+1,KK),BB1,BY2(KK),BY3(KK),LL,KFBI) 22002 CONTINUE DETADB1(KK)=PI2 ENDIF 12020 CONTINUE C FIP=FI+FSMM(J)+EPSM(J)*QCNTMI COSIJ=0. SINIJ=0. DO 2021KK=1,KMODXI DADU1(KK)=0. DADU2(KK)=0. IF(CSYM) GO TO 2021 DBDU1(KK)=0. DBDU2(KK)=0. 2021 CONTINUE DADU3=0. DADU4=0. DBDU3=0. DBDU4=0. DO 12021KK=1,KMODBI DADB1(KK)=0. DADB2(KK)=0. IF(CSYM) GO TO 12021 DBDB1(KK)=0. DBDB2(KK)=0. 12021 CONTINUE DADB3=0. DBDB3=0. DADB4=0. DBDB4=0. DO 2069M9=M9P,M9K DO 2068M8=M8P,M8K DO 2067M7=M7P,M7K DO 2066M6=M6P,M6K DO 2065M5=M5P,M5K DO 2060M4=M4P,M4K DO 2050M3=M3P,M3K DO 2040M2=M2P,M2K M1=MM-2*M2-3*M3-4*M4-M5-2*M6-3*M7-4*M8-M9 IF(KFXI.NE.0) THEN IF(KMODXI.EQ.2) M1=M1+M2 IF(KMODXI.EQ.3) M1=M1+2*M3 IF(KMODXI.EQ.4) M1=M1+3*M4 ENDIF IF(KFBI.NE.0) THEN IF(KMODBI.EQ.2) M1=M1+M6 IF(KMODBI.EQ.3) M1=M1+2*M7 IF(KMODBI.EQ.4) M1=M1+3*M8 ENDIF IF(M1.LT.M1P.OR.M1.GT.M1K) GO TO 2040 MS=M1+M2+M3+M4 IF(KFXI.NE.0) THEN IF(KMODXI.EQ.1) MS=MS-M1 IF(KMODXI.EQ.2) MS=MS-M2 IF(KMODXI.EQ.3) MS=MS-M3 IF(KMODXI.EQ.4) MS=MS-M4 ENDIF CHIS=0. S=1. IF(MOD(IABS(MS),2).EQ.1.AND..NOT.EPSL) S=-S IF(KK.EQ.KMODXI.AND.KFXI.NE.0.AND.EPSL) S=-S C DO 2030KK=1,KMODXI GO TO (2022,2023,2024,2025),KK 2022 FMP=-M1 MP=M1+5 GO TO 2026 2023 FMP=-M2 MP=M2+5 GO TO 2026 2024 FMP=-M3 MP=M3+5 GO TO 2026 2025 FMP=-M4 MP=M4+5 2026 IF(KK.EQ.KMODXI.AND.KFXI.NE.0) FMP=-FMP IF(EPSL) FMP=-FMP DSDU1(KK)=S*DBFDU1(MP,KK) DSDU2(KK)=S*DBFDU2(MP,KK) POM=BESP(MP,KK) S=S*POM IF(ABS(S).LT.DIF) GO TO 2040 IF(KK.LE.1) GO TO 2028 DO 2027LL=1,KK-1 DSDU1(LL)=DSDU1(LL)*POM 2027 DSDU2(LL)=DSDU2(LL)*POM 2028 DSDCHI(KK)=FMP CHIS=CHIS+FMP*CHI(KK) 2030 CONTINUE DO 12030KK=1,KMODBI GO TO (12022,12023,12024,12025),KK 12022 FMP=-M5 MP=IABS(M5)+1 GO TO 12026 12023 FMP=-M6 MP=IABS(M6)+1 GO TO 12026 12024 FMP=-M7 MP=IABS(M7)+1 GO TO 12026 12025 FMP=-M8 MP=IABS(M8)+1 12026 IF(KK.EQ.KMODBI.AND.KFBI.NE.0) FMP=-FMP IF(EPSL) FMP=-FMP DSDB1(KK)=S*DBFDB1(MP,KK) DSDB2(KK)=S*DBFDB2(MP,KK) IF(KK.EQ.KMODBI.AND.KFBI.NE.0) DSDB3(KK)=S*DBFDB3(MP,KK) POM=BESB(MP,KK) S=S*POM IF(ABS(S).LT.DIF) GO TO 2040 DO 12027LL=1,KK-1 DSDB1(LL)=DSDB1(LL)*POM IF(KK.EQ.KMODBI.AND.KFBI.NE.0) DSDB3(LL)=DSDB3(LL)*POM 12027 DSDB2(LL)=DSDB2(LL)*POM DO 22028LL=1,KMODXI DSDU1(LL)=DSDU1(LL)*POM 22028 DSDU2(LL)=DSDU2(LL)*POM 12028 DSDETA(KK)=FMP CHIS=CHIS+FMP*ETA(KK) IF(KK.NE.KMODBI.OR.KFBI.EQ.0) CHIS=CHIS-.25*FMP*PI2 12030 CONTINUE IF(KMODSI.LE.0) GO TO 2034 MP=IABS(M9) IF(MP.EQ.0) GO TO 2031 DSDAX=S*DSBDAX(MP) DSDAY=S*DSBDAY(MP) FMP=-SIGN(1.,FLOAT(M9)) IF(EPSL) FMP=-FMP CHIS=CHIS+KSI(MP)*FMP DSDKSI=FMP GO TO 12031 2031 DSDA0=S 12031 POM=SB(MP+1) S=S*POM DO 2032KK=1,KMODXI DSDU1(KK)=DSDU1(KK)*POM 2032 DSDU2(KK)=DSDU2(KK)*POM DO 2033KK=1,KMODBI DSDB1(KK)=DSDB1(KK)*POM IF(KK.EQ.KMODBI.AND.KFBI.NE.0) DSDB3(KK)=DSDB3(KK)*POM 2033 DSDB2(KK)=DSDB2(KK)*POM C 2034 POM=FIP+CHIS COSP=COS(POM) SINP=SIN(POM) COSPS=COSP*S SINPS=SINP*S COSIJ=COSIJ+COSPS SINIJ=SINIJ+SINPS DO 2035KK=1,KMODXI POMA=-DSDCHI(KK)*SINPS POMB= DSDCHI(KK)*COSPS IF(KK.NE.KMODXI.OR.KFXI.EQ.0) THEN DADU1(KK)=DADU1(KK)+POMA*DCHIDU1(KK)+COSP*DSDU1(KK) DADU2(KK)=DADU2(KK)+POMA*DCHIDU2(KK)+COSP*DSDU2(KK) ELSE DADU1(KK)=DADU1(KK)+COSP*DSDU1(KK) DADU2(KK)=DADU2(KK)+COSP*DSDU2(KK) DADU3 =DADU3 +COSP*DSDU3(KK) DADU4 =DADU4 +POMA*DCHIDU1(KK) ENDIF IF(CSYM) GO TO 2035 IF(KK.NE.KMODXI.OR.KFXI.EQ.0) THEN DBDU1(KK)=DBDU1(KK)+POMB*DCHIDU1(KK)+SINP*DSDU1(KK) DBDU2(KK)=DBDU2(KK)+POMB*DCHIDU2(KK)+SINP*DSDU2(KK) ELSE DBDU1(KK)=DBDU1(KK)+SINP*DSDU1(KK) DBDU2(KK)=DBDU2(KK)+SINP*DSDU2(KK) DBDU3 =DBDU3 +SINP*DSDU3(KK) DBDU4 =DBDU4 +POMB*DCHIDU1(KK) ENDIF 2035 CONTINUE DO 2036KK=1,KMODBI POMA=-DSDETA(KK)*SINPS POMB= DSDETA(KK)*COSPS IF(KK.NE.KMODBI.OR.KFBI.EQ.0) THEN DADB1(KK)=DADB1(KK)+POMA*DETADB1(KK)+COSP*DSDB1(KK) DADB2(KK)=DADB2(KK)+POMA*DETADB2(KK)+COSP*DSDB2(KK) ELSE DADB1(KK)=DADB1(KK)+COSP*DSDB1(KK) DADB2(KK)=DADB2(KK)+COSP*DSDB2(KK) DADB3 =DADB3 +COSP*DSDB3(KK) DADB4 =DADB4 +POMA*DETADB1(KK) ENDIF IF(CSYM) GO TO 2036 IF(KK.NE.KMODBI.OR.KFBI.EQ.0) THEN DBDB1(KK)=DBDB1(KK)+POMB*DETADB1(KK)+SINP*DSDB1(KK) DBDB2(KK)=DBDB2(KK)+POMB*DETADB2(KK)+SINP*DSDB2(KK) ELSE DBDB1(KK)=DBDB1(KK)+SINP*DSDB1(KK) DBDB2(KK)=DBDB2(KK)+SINP*DSDB2(KK) DBDB3 =DBDB3 +SINP*DSDB3(KK) DBDB4 =DBDB4 +POMB*DETADB1(KK) ENDIF 2036 CONTINUE IF(KMODSI.LE.0) GO TO 2040 IF(MP.EQ.0) GO TO 2038 POMA=-DSDKSI*SINPS POMB= DSDKSI*COSPS DADAX(MP)=DADAX(MP)+POMA*DKSIDAX(MP)+COSP*DSDAX DADAY(MP)=DADAY(MP)+POMA*DKSIDAY(MP)+COSP*DSDAY IF(CSYM) GO TO 2040 DBDAX(MP)=DBDAX(MP)+POMB*DKSIDAX(MP)+SINP*DSDAX DBDAY(MP)=DBDAY(MP)+POMB*DKSIDAY(MP)+SINP*DSDAY GO TO 2040 2038 DADA0=DADA0+DSDA0*COSP IF(CSYM) DBDA0=DBDA0+DSDA0*SINP 2040 CONTINUE 2050 CONTINUE 2060 CONTINUE 2065 CONTINUE 2066 CONTINUE 2067 CONTINUE 2068 CONTINUE 2069 CONTINUE C GO TO 2075 2070 SINIJ=SIN(FI) COSIJ=COS(FI) 2075 IF(IZO) GO TO 2076 EXPIJ=EXP(-HQJ*B1-KQJ*B2-LQJ*B3-HK2J*B4-HL2J*B5-KL2J*B6) COSIJ=COSIJ*EXPIJ SINIJ=SINIJ*EXPIJ C 2076 IF(HARM) GO TO 2080 DO 2077KK=1,KMODXI DADU1(KK)=EXPIJ*DADU1(KK) DADU2(KK)=EXPIJ*DADU2(KK) IF(KK.EQ.KMODXI.AND.KFXI.NE.0) THEN DADU3=EXPIJ*DADU3 DADU4=EXPIJ*DADU4 ENDIF IF(CSYM) GO TO 2077 DBDU1(KK)=EXPIJ*DBDU1(KK) DBDU2(KK)=EXPIJ*DBDU2(KK) IF(KK.EQ.KMODXI.AND.KFXI.NE.0) THEN DBDU3=EXPIJ*DBDU3 DBDU4=EXPIJ*DBDU4 ENDIF 2077 CONTINUE DO 12077KK=1,KMODBI DADB1(KK)=EXPIJ*DADB1(KK) DADB2(KK)=EXPIJ*DADB2(KK) IF(KK.EQ.KMODBI.AND.KFBI.NE.0) THEN DADB3=EXPIJ*DADB3 DADB4=EXPIJ*DADB4 ENDIF IF(CSYM) GO TO 12077 DBDB1(KK)=EXPIJ*DBDB1(KK) DBDB2(KK)=EXPIJ*DBDB2(KK) IF(KK.EQ.KMODBI.AND.KFBI.NE.0) THEN DBDB3=EXPIJ*DBDB3 DBDB4=EXPIJ*DBDB4 ENDIF 12077 CONTINUE DO 2078KK=1,KMODSI DADAX(KK)=EXPIJ*DADAX(KK) DADAY(KK)=EXPIJ*DADAY(KK) IF(CSYM) GO TO 2078 DBDAX(KK)=EXPIJ*DBDAX(KK) DBDAY(KK)=EXPIJ*DBDAY(KK) 2078 CONTINUE 2080 DADFII=DADFII+COSIJ IF(CSYM) GO TO 2090 DBDFII=DBDFII+SINIJ 2090 DADXI=DADXI-HJ*SINIJ DADYI=DADYI-KJ*SINIJ DADZI=DADZI-LJ*SINIJ IF(IZO) GO TO 2100 DADB1I=DADB1I-HQJ *COSIJ DADB2I=DADB2I-KQJ *COSIJ DADB3I=DADB3I-LQJ*COSIJ DADB4I=DADB4I-HK2J*COSIJ DADB5I=DADB5I-HL2J*COSIJ DADB6I=DADB6I-KL2J*COSIJ 2100 IF(CSYM) GO TO 2200 DBDXI=DBDXI+HJ*COSIJ DBDYI=DBDYI+KJ*COSIJ DBDZI=DBDZI+LJ*COSIJ IF(IZO) GO TO 2200 DBDB1I=DBDB1I-HQJ*SINIJ DBDB2I=DBDB2I-KQJ*SINIJ DBDB3I=DBDB3I-LQJ *SINIJ DBDB4I=DBDB4I-HK2J*SINIJ DBDB5I=DBDB5I-HL2J*SINIJ DBDB6I=DBDB6I-KL2J*SINIJ C 2200 IF(HARM) GO TO 5000 DO 3000KK=1,KMODXI DADUXX(KK)=DADUXX(KK)+HJ*DADU1(KK) DADUXY(KK)=DADUXY(KK)+KJ*DADU1(KK) DADUXZ(KK)=DADUXZ(KK)+LJ*DADU1(KK) IF(KK.NE.KMODXI.OR.KFXI.EQ.0) THEN DADUYX(KK)=DADUYX(KK)+HJ*DADU2(KK) DADUYY(KK)=DADUYY(KK)+KJ*DADU2(KK) DADUYZ(KK)=DADUYZ(KK)+LJ*DADU2(KK) ELSE DADUYX(KK)=DADUYX(KK)+DADU4 DADUYY(KK)=DADUYY(KK)+DADU2(KK) DADUYZ(KK)=DADUYZ(KK)+DADU3 ENDIF IF(CSYM) GO TO 3000 DBDUXX(KK)=DBDUXX(KK)+HJ*DBDU1(KK) DBDUXY(KK)=DBDUXY(KK)+KJ*DBDU1(KK) DBDUXZ(KK)=DBDUXZ(KK)+LJ*DBDU1(KK) IF(KK.NE.KMODXI.OR.KFXI.EQ.0) THEN DBDUYX(KK)=DBDUYX(KK)+HJ*DBDU2(KK) DBDUYY(KK)=DBDUYY(KK)+KJ*DBDU2(KK) DBDUYZ(KK)=DBDUYZ(KK)+LJ*DBDU2(KK) ELSE DBDUYX(KK)=DBDUYX(KK)+DBDU4 DBDUYY(KK)=DBDUYY(KK)+DBDU2(KK) DBDUYZ(KK)=DBDUYZ(KK)+DBDU3 ENDIF 3000 CONTINUE DO 4000KK=1,KMODBI DADBX1(KK)=DADBX1(KK)-HQJ*DADB1(KK) DADBX2(KK)=DADBX2(KK)-KQJ*DADB1(KK) DADBX3(KK)=DADBX3(KK)-LQJ*DADB1(KK) DADBX4(KK)=DADBX4(KK)-HK2J*DADB1(KK) DADBX5(KK)=DADBX5(KK)-HL2J*DADB1(KK) DADBX6(KK)=DADBX6(KK)-KL2J*DADB1(KK) IF(KK.NE.KMODBI.OR.KFBI.EQ.0) THEN DADBY1(KK)=DADBY1(KK)-HQJ*DADB2(KK) DADBY2(KK)=DADBY2(KK)-KQJ*DADB2(KK) DADBY3(KK)=DADBY3(KK)-LQJ*DADB2(KK) DADBY4(KK)=DADBY4(KK)-HK2J*DADB2(KK) DADBY5(KK)=DADBY5(KK)-HL2J*DADB2(KK) DADBY6(KK)=DADBY6(KK)-KL2J*DADB2(KK) ELSE DADBY1(KK)=DADBY1(KK)+DADB4 DADBY2(KK)=DADBY2(KK)+DADB2(KK) DADBY3(KK)=DADBY3(KK)+DADB3 ENDIF IF(CSYM) GO TO 4000 DBDBX1(KK)=DBDBX1(KK)-HQJ*DBDB1(KK) DBDBX2(KK)=DBDBX2(KK)-KQJ*DBDB1(KK) DBDBX3(KK)=DBDBX3(KK)-LQJ*DBDB1(KK) DBDBX4(KK)=DBDBX4(KK)-HK2J*DBDB1(KK) DBDBX5(KK)=DBDBX5(KK)-HL2J*DBDB1(KK) DBDBX6(KK)=DBDBX6(KK)-KL2J*DBDB1(KK) IF(KK.NE.KMODBI.OR.KFBI.EQ.0) THEN DBDBY1(KK)=DBDBY1(KK)-HQJ*DBDB2(KK) DBDBY2(KK)=DBDBY2(KK)-KQJ*DBDB2(KK) DBDBY3(KK)=DBDBY3(KK)-LQJ*DBDB2(KK) DBDBY4(KK)=DBDBY4(KK)-HK2J*DBDB2(KK) DBDBY5(KK)=DBDBY5(KK)-HL2J*DBDB2(KK) DBDBY6(KK)=DBDBY6(KK)-KL2J*DBDB2(KK) ELSE DBDBY1(KK)=DBDBY1(KK)+DBDB4 DBDBY2(KK)=DBDBY2(KK)+DBDB2(KK) DBDBY3(KK)=DBDBY3(KK)+DBDB3 ENDIF 4000 CONTINUE 5000 CONTINUE PH=FXI IF(PHFP.EQ.0..OR.PHF(I).EQ.0..OR.HARM) GO TO 5100 PH=FXI*EXP(PHFP*PHF(I)) 5100 IF(IZO) GO TO 5200 FTA=AII*PH FTACOS=FTA*DADFII DER(M )=PH*DADFII DER(M+4)=DADB1I*FTA DER(M+5)=DADB2I*FTA DER(M+6)=DADB3I*FTA DER(M+7)=DADB4I*FTA DER(M+8)=DADB5I*FTA DER(M+9)=DADB6I*FTA GO TO 5300 5200 TFI=EXP(-RHO*B1) FTA=PH*AII*TFI FTACOS=FTA*DADFII DADFII=TFI*DADFII DER(M )=PH*DADFII DER(M+4)=-RHO*FTACOS 5300 A=A+FTACOS AANOM=AANOM+FYI*FTACOS DER(M+1)=DADXI*FTA DER(M+2)=DADYI*FTA DER(M+3)=DADZI*FTA IF(HARM) GO TO 5400 IF(KMODSI.LE.0) GO TO 5311 DER(M+10)=DADA0*FTA LL=M+9 DO 5310KK=1,KMODSI LL=LL+2 DER(LL) =DADAX(KK)*FTA 5310 DER(LL+1)=DADAY(KK)*FTA 5311 LL=M+13 DO 5350KK=1,KMODXI LL=LL+6 DER(LL )=DADUXX(KK)*FTA DER(LL+1)=DADUXY(KK)*FTA DER(LL+2)=DADUXZ(KK)*FTA DER(LL+3)=DADUYX(KK)*FTA DER(LL+4)=DADUYY(KK)*FTA 5350 DER(LL+5)=DADUYZ(KK)*FTA LL=M+31 DO 5360KK=1,KMODBI LL=LL+12 DER(LL )=DADBX1(KK)*FTA DER(LL+1 )=DADBX2(KK)*FTA DER(LL+2 )=DADBX3(KK)*FTA DER(LL+3 )=DADBX4(KK)*FTA DER(LL+4 )=DADBX5(KK)*FTA DER(LL+5 )=DADBX6(KK)*FTA DER(LL+6 )=DADBY1(KK)*FTA DER(LL+7 )=DADBY2(KK)*FTA DER(LL+8 )=DADBY3(KK)*FTA DER(LL+9 )=DADBY4(KK)*FTA DER(LL+10)=DADBY5(KK)*FTA 5360 DER(LL+11)=DADBY6(KK)*FTA DER(M+91)=PHFP*FTACOS 5400 IF(CSYM) GO TO 6000 FTASIN=FTA*DBDFII B=B+FTASIN BANOM=BANOM+FYI*FTASIN DC(M+1)=DBDXI*FTA DC(M+2)=DBDYI*FTA DC(M+3)=DBDZI*FTA IF(IZO) GO TO 5600 DC(M )=PH*DBDFII DC(M+4 )=DBDB1I*FTA DC(M+5 )=DBDB2I*FTA DC(M+6 )=DBDB3I*FTA DC(M+7 )=DBDB4I*FTA DC(M+8 )=DBDB5I*FTA DC(M+9 )=DBDB6I*FTA GO TO 5700 5600 DBDFII=TFI*DBDFII DC(M )=PH*DBDFII DC(M+4 )=-RHO*FTASIN 5700 IF(HARM) GO TO 6000 DC(M+10)=DBDA0*FTA LL=M+9 DO 5710KK=1,KMODSI LL=LL+2 DC(LL) =DBDAX(KK)*FTA 5710 DC(LL+1)=DBDAY(KK)*FTA LL=M+13 DO 5750KK=1,KMODXI LL=LL+6 DC(LL )=DBDUXX(KK)*FTA DC(LL+1)=DBDUXY(KK)*FTA DC(LL+2)=DBDUXZ(KK)*FTA DC(LL+3)=DBDUYX(KK)*FTA DC(LL+4)=DBDUYY(KK)*FTA 5750 DC(LL+5)=DBDUYZ(KK)*FTA LL=M+31 DO 5460KK=1,KMODBI LL=LL+12 DC(LL )=DBDBX1(KK)*FTA DC(LL+1 )=DBDBX2(KK)*FTA DC(LL+2 )=DBDBX3(KK)*FTA DC(LL+3 )=DBDBX4(KK)*FTA DC(LL+4 )=DBDBX5(KK)*FTA DC(LL+5 )=DBDBX6(KK)*FTA DC(LL+6 )=DBDBY1(KK)*FTA DC(LL+7 )=DBDBY2(KK)*FTA DC(LL+8 )=DBDBY3(KK)*FTA DC(LL+9 )=DBDBY4(KK)*FTA DC(LL+10)=DBDBY5(KK)*FTA 5460 DC(LL+11)=DBDBY6(KK)*FTA DC(M+91)=PHFP*FTASIN 6000 M=M+92 A=A-BANOM B=B+AANOM A=A*CENTR B=B*CENTR SQRTAB=SQRT(A**2+B**2) DER(IQ)=SQRTAB YC=SC(IQ)*SQRTAB IF(SQRTAB.EQ.0.) GO TO 6100 SQRTAB=SC(IQ)/SQRTAB 6100 C1=A*SQRTAB*CENTR C2=B*SQRTAB*CENTR M=18 IF(CSYM) GO TO 9500 DO 9200I=1,NA KMODXI=KMODX(I) KMODBI=KMODB(I) KMODSI=KMODS(I) IF(KMODXI.LT.0) KMODXI=KMODXM(-KMODXI) J=ISF(I) FYI=FY(J) CA=C1+FYI*C2 CB=C2-FYI*C1 J=ITF(I)*5 DO 9100N=1,J M=M+1 9100 DER(M)=CA*DER(M)+CB*DC(M) IF(J.NE.10) M=M+5 IF(KMODSI.GT.0) GO TO 9120 M=M+9 GO TO 9130 9120 M=M+1 DER(M)=CA*DER(M)+CB*DC(M) KK=KMODSI*2 DO 9125N=1,KK M=M+1 9125 DER(M)=CA*DER(M)+CB*DC(M) M=M+8-KK 9130 IF(KMODXI.GT.0) GO TO 9140 M=M+24 GO TO 9160 9140 KK=KMODXI*6 DO 9150N=1,KK M=M+1 9150 DER(M)=CA*DER(M)+CB*DC(M) M=M+24-KK 9160 IF(KMODBI.GT.0) GO TO 9170 M=M+49 GO TO 9190 9170 KK=KMODBI*12 DO 9180N=1,KK M=M+1 9180 DER(M)=CA*DER(M)+CB*DC(M) M=M+49-KK 9190 DER(M)=CA*DER(M)+CB*DC(M) 9200 CONTINUE GO TO 9900 9500 DO 9300I=1,NA KMODXI=KMODX(I) KMODBI=KMODB(I) KMODSI=KMODS(I) IF(KMODXI.LT.0) KMODXI=KMODXM(-KMODXI) J=ISF(I) CA=C1+FY(J)*C2 J=ITF(I)*5 DO 9400N=1,J M=M+1 9400 DER(M)=CA*DER(M) IF(J.NE.10) M=M+5 IF(KMODSI.GT.0) GO TO 9420 M=M+9 GO TO 9430 9420 M=M+1 DER(M)=CA*DER(M) KK=KMODSI*2 DO 9425N=1,KK M=M+1 9425 DER(M)=CA*DER(M) M=M+8-KK 9430 IF(KMODXI.GT.0) GO TO 9440 M=M+24 GO TO 9460 9440 KK=KMODXI*6 DO 9450N=1,KK M=M+1 9450 DER(M)=CA*DER(M) M=M+24-KK 9460 IF(KMODBI.GT.0) GO TO 9470 M=M+49 GO TO 9490 9470 KK=KMODBI*12 DO 9480N=1,KK M=M+1 9480 DER(M)=CA*DER(M) M=M+49-KK 9490 DER(M)=CA*DER(M) 9300 CONTINUE C 9900 IF(IEXT.LE.0) GO TO 9999 CALL PHILIP(YC,EXTINC,EXTINK,DER(7),SINTHL,SC(IQ)) EXTKOR=SQRT(EXTINC) YC=YC*EXTKOR DER(IQ)=DER(IQ)*EXTKOR FI=EXTKOR*EXTINK/EXTINC DO 9910I=19,NP 9910 DER(I)=DER(I)*FI 9999 RETURN END SUBROUTINE PHILIP(F,EXTINC,EXTINK,DFDC,SINTHL,SC) DIMENSION DFDC(12) COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10), 1 BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ, 2 EC(12) YY(A,B,C,FX)=1./(1.+C*FX+A*FX**2/(1.+B*FX)) ZZ(A,B,C,FX,FY)=-0.5*FY*(C+A*FX*(2.+B*FX)/(1.+B*FX)**2) DO 30I=1,12 30 DFDC(I)=0. IF(F.LE.0.0) RETURN FSQR=(F/SC)**2 TBAR=EFPIP(7)/VOLUME GAMMA=1.5E-4*FLAME/VOLUME HROOT=EC(1) IF(ITYPEX.EQ.1) HROOT=EC(7) IF(IEXT.EQ.1) GO TO 4 N=0 M=6 IF(IEXT.NE.3) GO TO 3 N=3 M=0 3 HROOT=EFPIP(N+1)**2*EC(M+1)+EFPIP(N+2)**2*EC(M+2)+EFPIP(N+3)**2 1 *EC(M+3)+2.*(EFPIP(N+1)*EFPIP(N+2)*EC(M+4)+EFPIP(N+1)* 2 EFPIP(N+3)*EC(M+5)+EFPIP(N+2)*EFPIP(N+3)*EC(M+6)) HROOT=SQRT(HROOT) IF(IEXT.EQ.3.OR.IMOSAIC.EQ.1) HROOT=1./HROOT 4 SINT=SINTHL*FLAME TH=ASIN(SINT)/TORAD COST=SQRT(1.-SINT**2) IF(IAB.EQ.0) GO TO 41 I=TH*.2 XI=TH*.2-FLOAT(I) I=I+1 CP=.6666667*(DADMI(I)+(DADMI(I+1)-DADMI(I))*XI) TBAR=TBAR*CP GAMMA=GAMMA*CP 41 SIN2T=2.*SINT*COST COS2T=COST**2-SINT**2 AP=0.20+0.45*COS2T BP=0.22-0.12*(0.5-COS2T)**2 CP=2. IF(IAB.EQ.0) GO TO 42 I=SINT*10. XI=SINT*10.-FLOAT(I) I=I+1 AG=AA(I)+(AA(I+1)-AA(I))*XI BG=BA(I)+(BA(I+1)-BA(I))*XI AL=AG BL=BG GO TO 43 42 AG=0.58+0.48*COS2T+0.24*COS2T**2 BG=0.020-0.025*COS2T AL=0.025+0.285*COS2T BL=0.15-0.2*(0.75-COS2T)**2 IF(COS2T.LT.0.) BL=-0.45*COS2T 43 IF(ITYPEX.NE.2) GO TO 130 AS=AP BS=BP CS=2. GO TO 140 130 AS=AG BS=BG CS=2.12 IF(IDISTR.EQ.1) GO TO 140 AS=AL BS=BL CS=2. 140 IF(ITYPEX.NE.2) GO TO 100 C HERE FOR TYPE 2 RHO=HROOT BETA=RHO DBRHO=1. GO TO 120 100 IF(ITYPEX.NE.1)GO TO 110 C HERE FOR TYPE 1 G=HROOT BETA=G/SIN2T DBG=1./SIN2T RHO=0. DBRHO=0. GO TO 120 C HERE FOR GENERAL CASE 110 RHO=EC(1) G=EC(7) IF(IEXT.EQ.1) GO TO 150 G=HROOT 150 BETA=RHO/SQRT(1.+(RHO*SIN2T/G)**2) DBRHO=(BETA/RHO)**3 DBG=(BETA/G)**3*SIN2T**2 IF(IDISTR.EQ.1)GO TO 120 BETA=RHO/(1.+RHO*SIN2T/G) DBRHO=(BETA/RHO)**2 DBG=(BETA/G)**2*SIN2T 120 CONTINUE X2=TBAR*BETA*FSQR X1=0. C INTRODUCE PRIMARY EXTINCTION HERE FOR GENERAL CASE IF(ITYPEX.EQ.3) X1=GAMMA*FSQR*RHO**2 XP=X1 YP=YY(AP,BP,CP,XP) ZP=ZZ(AP,BP,CP,XP,YP) YP=SQRT(YP) IF(ITYPEX.NE.3) ZP=0. XS=X2*YP YS=YY(AS,BS,CS,XS) ZS=ZZ(AS,BS,CS,XS,YS) YS=SQRT(YS) C XRAY CASE TREATED HERE POL=COS2T**2 FACT=1./(C2THMQ+POL) XPP=X1*POL YPP=YY(AP,BP,CP,XPP) ZPP=ZZ(AP,BP,CP,XPP,YPP) IF(ITYPEX.NE.3) ZPP=0. YPP=SQRT(YPP) XPS=YPP*X2*POL YPS=YY(AS,BS,CS,XPS) ZPS=ZZ(AS,BS,CS,XPS,YPS) YPS=SQRT(YPS) YPYS=C2THMQ*YP*YS YPPYPS=POL*YPP*YPS EXTINC=FACT*(YPYS+YPPYPS) OM1=FACT*(YPYS*(1.+XP*ZP+XS*ZS*(1.+XP*ZP))+YPPYPS* 1 (1.+XPP*ZPP+XPS*ZPS*(1.+XPP*ZPP))) OM2=FACT*(YPYS*ZS*YP+POL*YPPYPS*ZPS*YPP) OM3=FACT*(YPYS*(2.*GAMMA*RHO*ZP+ZS*YP*(TBAR*DBRHO +2.*X2*GAMMA* 1 RHO*ZP))+POL*YPPYPS*(2.*GAMMA*RHO*ZPP+ZPS*YPP*(TBAR*DBRHO 2 +2.*POL*X2*GAMMA*RHO*ZPP))) EXTINK=OM1 TUMS=0.5*(SC*FSQR)**2/(F*SQRT(EXTINC)) TUMS1=TUMS*OM3 TUMS2=TUMS*TBAR*OM2*DBG DFDC(1)=TUMS1 IF(IEXT.NE.1) GO TO 7 DFDC(7)=TUMS2 GO TO 70 7 DIVID=0.5/HROOT IF(IEXT.EQ.3.OR.IMOSAIC.EQ.1) DIVID=-0.5*HROOT**3 IF(ITYPEX.EQ.2) TUMS2=TUMS1 DFDC(M+1)=TUMS2*EFPIP(N+1)**2*DIVID DFDC(M+2)=TUMS2*EFPIP(N+2)**2*DIVID DFDC(M+3)=TUMS2*EFPIP(N+3)**2*DIVID DFDC(M+4)=2.*TUMS2*EFPIP(N+1)*EFPIP(N+2)*DIVID DFDC(M+5)=2.*TUMS2*EFPIP(N+1)*EFPIP(N+3)*DIVID DFDC(M+6)=2.*TUMS2*EFPIP(N+2)*EFPIP(N+3)*DIVID 70 RETURN END SUBROUTINE SUMA COMMON/REFL/ YO,YC,A,B,DY,WDYQ,NULOVA,SINTHL,RHO,WT,DER(9708), 1 DC(9708),DYP,IQ COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10), 1 UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10), 2 PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10), 3 SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100), 4 NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10), 5 UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10) COMMON NPS(9),NAM(9),NPSS(10),NAMM,NP,NZ,TLUM,Q,KI(9708),PS(9708), 1 AM(75000) LOGICAL NULOVA DIMENSION IPSA(9) DO 500IP=1,9 500 IPSA(IP)=NPSS(IP) DO 700IP=1,NPMK N=KI(IP) IF(N.EQ.0) GO TO 700 IPS=IPSA(N)+1 DC(IPS)=DER(IP) IPSA(N)=IPS 700 CONTINUE IM=0 IPK=0 DO 4000N=1,9 NPSN=NPS(N) IF(NPSN.LE.0) GO TO 4000 IPZ=IPK+1 IPK=IPZ+NPSN-1 DO 3000IP=IPZ,IPK DERI=DC(IP)*WT IF(DERI.NE.0.0) GO TO 1000 IM=IM+IP-IPZ+1 GO TO 3000 1000 DO 2000JP=IPZ,IP IM=IM+1 DERJ=DC(JP) IF(DERJ.EQ.0.) GO TO 2000 AM(IM)=AM(IM)+DERI*DERJ 2000 CONTINUE PS(IP)=PS(IP)+DERI*DYP 3000 CONTINUE 4000 CONTINUE RETURN END SUBROUTINE INVER C COMMON/KLICE/ NA,NS,NF,CSYM,VYH,NC,NOWR,IWQ,YOMIN,YOMAX,BLKOEF, 1 CENTR,NCYKL,SNLMN,SNLMX,KIM,DIF,KIW,LITE,IFSQ COMMON/UNITS/ W,M40,M50 COMMON/REFL/ YO,YC,A,B,DY,WDYQ,NULOVA,SINTHL,RHO,WT,DER(9708), 1 DC(9708),DYP,IQ COMMON/CM4/ ATOM(100),AI(100),ISF(100),ITF(100),X(3,100), 1 BETA(6,100),SC(6),KMODX(100),UX(3,4,100),UY(3,4,100) 2 ,PHF(100),QCNTM(100),KMODS(100),A0(100),AX(4,100), 3 AY(4,100),KMODB(100),BX(6,4,100),BY(6,4,100), 4 KFX(100),KFB(100) COMMON/CSM4/ SAI(100),SX(3,100),SB(6,100),SUX(3,4,100), 1 SUY(3,4,100),SPHF(100),SA0(100),SAX(4,100), 2 SAY(4,100),SBX(6,4,100),SBY(6,4,100) COMMON/LABEL/ L1,L2(3),L5,L7,L8,L9,L10 COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10), 1 BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ, 2 EC(12) COMMON/VAZBY/ NVX,NX(10),NVB,NB(10),RMX(3,3,10),RMB(6,6,10), 1 NVAI,NAI(2,10),NAIXB(10),SUMAI(10) COMMON/MODUL/ QUI(3),QUIR(3) COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10), 1 UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10), 2 PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10), 3 SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100), 4 NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10), 5 UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10) COMMON/STALY/ PI2 COMMON NPS(9),NAM(9),NPSS(10),NAMM,NP,NZ,TLUM,Q,KI(9708),PS(9708), 1 AM(75000) EQUIVALENCE (NPSS(10),NPSM) DIMENSION ISING(9),ECP(12),SCP(6) LOGICAL CSYM INTEGER W CHARACTER*1 CU,CV,CW CHARACTER*8 ATOM,AT,MOLNAME CHARACTER*12 L1,L2,L5,L7,L8,L9,L10,LP DATA AT/' '/,CU,CV,CW/'U','V','W'/ IPS=1 IAM=1 Q=Q/FLOAT(NZ-NPSM-1) DO 3000I=1,9 N=NPS(I) IF(N.LE.0) GO TO 3000 C CALL SMI(AM(IAM),N,ISING(I),I) CALL SKOVEJ(AM(IAM),Q,I,KI) IF(ISING(I).EQ.0) GO TO 1500 C KK=IPS+N-1 DO 1200K=IPS,KK 1200 DC(K)=0. WRITE(W,101) I GO TO 2800 1500 CALL NASOB(AM(IAM),PS(IPS),DC(IPS),N) 2800 IAM=IAM+NAM(I) IPS=IPS+N 3000 CONTINUE Q=SQRT(Q) DO 3100I=1,NPMK PS(I)=0. 3100 DER(I)=0. IAM=0 IPS=0 DO 3500I=1,9 N=NPS(I) IF(N.EQ.0) GO TO 3500 IF(ISING(I).EQ.0) GO TO 3200 IAM=IAM+NAM(I) IPS=IPS+N GO TO 3500 3200 K=0 C DO 3300IP=1,NPMK IF(KI(IP).NE.I) GO TO 3300 K=K+1 IAM=IAM+K IPS=IPS+1 PS(IP)=DC(IPS) DER(IP)=SQRT(AM(IAM))*Q 3300 CONTINUE 3500 CONTINUE IF(NVX.LE.0) GO TO 3525 DO 3521I=1,NVX K=NX(I)*92-72 3521 CALL ZMZPET(PS(K),DER(K),RMX(1,1,I),3) 3525 IF(NVB.LE.0) GO TO 3529 DO 3526I=1,NVB K=NB(I)*92-69 3526 CALL ZMZPET(PS(K),DER(K),RMB(1,1,I),6) 3529 ISPASA=0 TLUMP=TLUM DO 3530I=1,6 3530 SCP(I)=SC(I) DO 3501I=1,12 3501 ECP(I)=EC(I) 3503 CALL ZMENA(L1,POM,POM,1,IT,AT,IP) IT=1 CALL PIS(IT) DO 3504I=1,6 ENCODE(12,100,LP) I 3504 CALL ZMENA(LP,SC(I),POM,2,IT,AT,IP) IF(IEXT.LE.0) GO TO 3690 IF(ITYPEX.EQ.1) GO TO 3600 IF(IEXT.EQ.3) GO TO 3550 CALL ZMENA(L7,EC(1),POM,2,IT,AT,IP) IF(EC(1).LT.0.) GO TO 3680 IF(ITYPEX.EQ.3) GO TO 3600 GO TO 3690 3550 DO 3560J=1,6 CALL INDEX(J,K,L) ENCODE(12,110,LP) K,L 3560 CALL ZMENA(LP,EC(J),POM,2,IT,AT,IP) IF(EC(1).LE.0.) GO TO 3680 IF(EC(1)*EC(2)-EC(4)**2.LE.0.) GO TO 3680 IF(EC(1)*EC(2)*EC(3)+2.*EC(4)*EC(5)*EC(6)-EC(1)*EC(6)**2 1 -EC(2)*EC(5)**2-EC(3)*EC(4)**2.LE.0.) GO TO 3680 GO TO 3690 3600 IP=12 IF(IEXT.GT.1) GO TO 3650 CALL ZMENA(L8,EC(7),POM,2,IT,AT,IP) IF(EC(7).LE.0.) GO TO 3680 GO TO 3690 3650 DO 3660J=1,6 CALL INDEX(J,K,L) ENCODE(12,111,LP) K,L 3660 CALL ZMENA(LP,EC(J+6),POM,2,IT,AT,IP) IF(EC(7).LE.0.) GO TO 3680 IF(EC(7)*EC(8)-EC(10)**2.LE.0.) GO TO 3680 IF(EC(7)*EC(8)*EC(9)+2.*EC(10)*EC(11)*EC(12)-EC(7)*EC(12)**2 1 -EC(8)*EC(11)**2-EC(9)*EC(10)**2.LE.0.) GO TO 3680 GO TO 3690 3680 ISPASA=ISPASA+1 DO 13680I=1,6 13680 SC(I)=SCP(I) DO 3681I=1,12 3681 EC(I)=ECP(I) TLUM=TLUM*.5 IF(ISPASA.GT.5) STOP WRITE(W,109) TLUMP,TLUM GO TO 3503 C C .... CALCULATION OF CHANGES OF ATOMIC PARAMETERS C 3690 IP=18 DO 4000I=1,NA DO 3691K=1,NVAI IF(I.EQ.NAI(2,K)) GO TO 3692 3691 CONTINUE GO TO 3699 3692 N=NAI(1,K) IF(NAIXB(K).LT.0) GO TO 13692 AI(I)=SUMAI(K)-AI(N) GO TO 33692 13692 A0(I)=SUMAI(K)-A0(N) DO 23692KK=1,KMODS(I) AX(KK,I)=-AX(KK,N) 23692 AY(KK,I)=-AY(KK,N) 33692 IF(IABS(NAIXB(K)).EQ.3) GO TO 3699 DO 3693KK=1,3 3693 X(KK,I)=X(KK,N) IF(IABS(NAIXB(K)).EQ.2) GO TO 3699 DO 3694KK=1,6 3694 BETA(KK,I)=BETA(KK,N) DO 3696KK=1,KMODX(I) DO 3695L=1,3 UX(L,KK,I)=UX(L,KK,N) 3695 UY(L,KK,I)=UY(L,KK,N) 3696 CONTINUE DO 3698KK=1,KMODB(I) DO 3697L=1,6 BX(L,KK,I)=BX(L,KK,N) 3697 BY(L,KK,I)=BY(L,KK,N) 3698 CONTINUE PHF(I)=PHF(N) C 3699 CALL PIS(IT) WRITE(W,102) ATOM(I) CALL ZMENA(L1,AI(I),SAI(I),2,IT,ATOM(I),IP) CALL ZMENA(L2(1),X(1,I),SX(1,I),2,IT,ATOM(I),IP) CALL ZMENA(L2(2),X(2,I),SX(2,I),2,IT,ATOM(I),IP) CALL ZMENA(L2(3),X(3,I),SX(3,I),2,IT,ATOM(I),IP) IF(ITF(I).NE.1) GO TO 3700 CALL ZMENA(L5,BETA(1,I),SB(1,I),2,IT,ATOM(I),IP) IP=IP+5 GO TO 3800 3700 DO 3710J=1,6 CALL INDEX(J,K,L) IF(LITE.EQ.0) GO TO 3705 ENCODE(12,112,LP) K,L GO TO 3710 3705 ENCODE(12,117,LP) K,L 3710 CALL ZMENA(LP,BETA(J,I),SB(J,I),2,IT,ATOM(I),IP) C 3800 IF(KMODS(I).LE.0) GO TO 3808 CALL ZMENA(L10,A0(I),SA0(I),2,IT,ATOM(I),IP) DO 3802J=1,4 IF(J.GT.KMODS(I)) GO TO 3801 ENCODE(12,118,LP) J CALL ZMENA(LP,AX(J,I),SAX(J,I),2,IT,ATOM(I),IP) ENCODE(12,119,LP) J CALL ZMENA(LP,AY(J,I),SAY(J,I),2,IT,ATOM(I),IP) GO TO 3802 3801 IP=IP+2 3802 CONTINUE GO TO 3809 3808 IP=IP+9 3809 IF(KMODX(I).LE.0.AND.KMODB(I).LE.0) GO TO 3995 DO 3900J=1,4 IF(J.GT.KMODX(I)) GO TO 3850 DO 3810K=1,3 ENCODE(12,113,LP) CU,K,J 3810 CALL ZMENA(LP,UX(K,J,I),SUX(K,J,I),2,IT,ATOM(I),IP) DO 3820K=1,3 ENCODE(12,114,LP) CU,K,J 3820 CALL ZMENA(LP,UY(K,J,I),SUY(K,J,I),2,IT,ATOM(I),IP) GO TO 3900 3850 IP=IP+6 3900 CONTINUE DO 3990J=1,4 IF(J.GT.KMODB(I)) GO TO 3950 DO 3910K=1,6 ENCODE(12,113,LP) 'T',K,J 3910 CALL ZMENA(LP,BX(K,J,I),SBX(K,J,I),2,IT,ATOM(I),IP) DO 3920K=1,6 ENCODE(12,114,LP) 'T',K,J 3920 CALL ZMENA(LP,BY(K,J,I),SBY(K,J,I),2,IT,ATOM(I),IP) GO TO 3990 3950 IP=IP+12 3990 CONTINUE CALL ZMENA(L9,PHF(I),SPHF(I),2,IT,ATOM(I),IP) GO TO 4000 3995 IP=IP+73 4000 CONTINUE CALL RESETX(X) CALL RESETB(BETA) CALL RESETA(A0,AX,AY) CALL RESETD(UX,UY) CALL RESETC(BX,BY) IF(NMOL.EQ.0) GO TO 4200 DO 4100I=1,NMOL CALL PIS(IT) WRITE(W,102) MOLNAME(I) IF(KMODXM(I).LE.0) GO TO 4090 DO 4040J=1,4 IF(J.GT.KMODXM(I)) GO TO 4030 DO 4010K=1,3 ENCODE(12,113,LP) CV,K,J 4010 CALL ZMENA(LP,UTX(K,J,I),SUTX(K,J,I),2,IT,MOLNAME(I),IP) DO 4020K=1,3 ENCODE(12,114,LP) CV,K,J 4020 CALL ZMENA(LP,UTY(K,J,I),SUTY(K,J,I),2,IT,MOLNAME(I),IP) GO TO 4040 4030 IP=IP+6 4040 CONTINUE DO 4080J=1,4 IF(J.GT.KMODXM(I)) GO TO 4070 DO 4050K=1,3 ENCODE(12,113,LP) CW,K,J 4050 CALL ZMENA(LP,URX(K,J,I),SURX(K,J,I),2,IT,MOLNAME(I),IP) DO 4060K=1,3 ENCODE(12,114,LP) CW,K,J 4060 CALL ZMENA(LP,URY(K,J,I),SURY(K,J,I),2,IT,MOLNAME(I),IP) GO TO 4080 4070 IP=IP+6 4080 CONTINUE CALL ZMENA(L9,PHFM(I),SPHFM(I),2,IT,MOLNAME(I),IP) GO TO 4100 4090 IP=IP+49 4100 CONTINUE CALL RESETM(UTX,UTY,URX,URY) 4200 CALL ZMENA(L1,POM,POM,3,IT,AT,IP) DO 4500I=1,NA IF(BETA(1,I).GT.0.) GO TO 4510 WRITE(W,107) ATOM(I) GO TO 4500 4510 IF(ITF(I).EQ.1) GO TO 4500 IF(BETA(1,I)*BETA(2,I)-BETA(4,I)**2.GT.0.) GO TO 4520 WRITE(W,107) ATOM(I) GO TO 4500 4520 IF(BETA(1,I)*BETA(2,I)*BETA(3,I)+2.*BETA(4,I)*BETA(5,I)*BETA(6,I) 1 -BETA(3,I)*BETA(4,I)**2-BETA(1,I)*BETA(6,I)**2-BETA(2,I)* 2 BETA(5,I)**2.GT.0.) GO TO 4500 WRITE(W,107) ATOM(I) 4500 CONTINUE REWIND M40 POM=0.0 WRITE(M40,103) NA,NMOL WRITE(M40,108) SC,(KI(I),I=1,6) WRITE(M40,104) POM WRITE(M40,108)(EC(I),I=1,6),(KI(I),I=7,12), 1 (EC(I),I=7,12),(KI(I),I=13,18) K=18 DO 5000I=1,NA KK=K+1 K=K+10 J=IABS(KMODX(I))+10*KMODB(I) IF(KMODX(I).LT.0) J=-J WRITE(M40,105) ATOM(I),ISF(I),KFB(I),KFX(I),ITF(I),J,KMODS(I), 1 AI(I),(X(N,I),N=1,3),(BETA(N,I),N=1,6), 2 (KI(N),N=KK,K) KK=K+1 K=K+1 IF(KMODS(I).GT.0) WRITE(M40,104) A0(I),KI(K) DO 4800J=1,4 KK=K+1 K=K+2 IF(J.GT.KMODS(I)) GO TO 4800 WRITE(M40,120) AX(J,I),AY(J,I),KI(KK),KI(K) 4800 CONTINUE DO 4900J=1,4 KK=K+1 K=K+6 IF(J.GT.KMODX(I)) GO TO 4900 WRITE(M40,108)(UX(N,J,I),N=1,3),(UY(N,J,I),N=1,3),(KI(N),N=KK,K) 4900 CONTINUE DO 4950J=1,4 KK=K+1 K=K+12 IF(J.GT.KMODB(I)) GO TO 4950 WRITE(M40,108)(BX(N,J,I),N=1,6),(KI(N),N=KK,KK+5), 1 (BY(N,J,I),N=1,6),(KI(N),N=KK+6,K) 4950 CONTINUE KK=K+1 K=K+1 IF(KMODX(I)+KMODS(I)+KMODB(I).GT.0) WRITE(M40,104) PHF(I),KI(KK) 5000 CONTINUE IF(NMOL.EQ.0) GO TO 5200 DO 5100I=1,NMOL CALL CENTER(I) QCNTMM=-PI2*((QUI(1)+QUIR(1))*XM(1,I)+(QUI(2)+QUIR(2))*XM(2,I)+ 1 (QUI(3)+QUIR(3))*XM(3,I)) DO 5003N=1,NA IF(KMODX(N).NE.-I) GO TO 5003 QCNTM(N)=QCNTMM 5003 CONTINUE 5005 WRITE(M40,116) MOLNAME(I),NPG(I),KMODXM(I),(XM(N,I),N=1,3) IF(NPG(I).GT.1) WRITE(M40,103)(KPG(J,I),J=2,NPG(I)) DO 5010J=1,4 KK=K+1 K=K+6 IF(J.GT.KMODXM(I)) GO TO 5010 WRITE(M40,108)(UTX(N,J,I),N=1,3),(UTY(N,J,I),N=1,3),(KI(N),N=KK,K) 5010 CONTINUE DO 5020J=1,4 KK=K+1 K=K+6 IF(J.GT.KMODXM(I)) GO TO 5020 WRITE(M40,108)(URX(N,J,I),N=1,3),(URY(N,J,I),N=1,3),(KI(N),N=KK,K) 5020 CONTINUE KK=K+1 K=K+1 IF(KMODXM(I).GT.0) WRITE(M40,104) PHFM(I),KI(KK) 5100 CONTINUE 5200 DO 6000I=1,NA WRITE(M40,106) ATOM(I),SAI(I),(SX(N,I),N=1,3),(SB(N,I),N=1,6) K=KMODS(I) IF(K.LE.0) GO TO 5400 WRITE(M40,104) SA0(I) DO 5300J=1,K 5300 WRITE(M40,120) SAX(J,I),SAY(J,I) 5400 DO 5500J=1,KMODX(I) 5500 WRITE(M40,108)(SUX(N,J,I),N=1,3),(SUY(N,J,I),N=1,3) DO 5600J=1,KMODB(I) WRITE(M40,108)(SBX(N,J,I),N=1,6) 5600 WRITE(M40,108)(SBY(N,J,I),N=1,6) IF(KMODX(I)+KMODS(I)+KMODB(I).GT.0) WRITE(M40,108) SPHF(I) 6000 CONTINUE IF(NMOL.EQ.0) GO TO 7000 DO 6100I=1,NMOL K=KMODXM(I) WRITE(M40,106) MOLNAME(I) IF(K.LE.0) GO TO 6100 DO 6010J=1,K 6010 WRITE(M40,108)(SUTX(N,J,I),N=1,3),(SUTY(N,J,I),N=1,3) DO 6020J=1,K 6020 WRITE(M40,108)(SURX(N,J,I),N=1,3),(SURY(N,J,I),N=1,3) WRITE(M40,108) SPHFM(I) 6100 CONTINUE 7000 IF(ISPASA.NE.0) TLUM=TLUMP 100 FORMAT('SCALE-',I1) 101 FORMAT(' MATRIX',I2,'.BLOKU IS SINGULAR') 102 FORMAT('0',A8) 103 FORMAT(12I5) 104 FORMAT(F9.6,51X,I1) 105 FORMAT(A8,I3,3I1,I3,I1,4F9.6/6F9.6,6X,10I1) 106 FORMAT(A8,10X,4F9.6/6F9.6) 107 FORMAT('TTEPMERATURE FACTOR ATOMU : ',A8,' IS NOT POSITIVE DEFINIT 1E') 108 FORMAT(6F9.6,6X,6I1) 109 FORMAT('0*** BEWARE ||| ***'/'0DUMPING FACTOR WAS CHANGED FROM' 1,F10.4,' TO',F10.4/' IN ORDER TO SAVE POSITIVE DEFINITY OF EXTINCT 2IONAL PARAMETERS') 110 FORMAT('RHO(',I1,',',I1,') ') 111 FORMAT('G(',I1,',',I1,') ') 112 FORMAT('BETA(',I1,',',I1,') ') 113 FORMAT(A1,I1,'*SIN(',I1,'*X4)') 114 FORMAT(A1,I1,'*COS(',I1,'*X4)') 116 FORMAT(A8,3X,2I3,10X,3F9.6) 117 FORMAT('U(',I1,',',I1,') ') 118 FORMAT('A*SIN(',I1,'*X4) ') 119 FORMAT('A*COS(',I1,'*X4) ') 120 FORMAT(2F9.6,42X,2I1) RETURN END SUBROUTINE PIS(IT) COMMON/UNITS/ W,M40,M50 INTEGER W IF(MOD(IT,50).NE.1) GO TO 9999 CALL NEWPG WRITE(W,100) 100 FORMAT('0PRINT OF CHANGES OF PARAMETERS :'/'0PARAMETER',9X,' OLD 1' 10X,'CHANGE',7X,'NEW',10X,'E.S.D.',1X,'CHANGE/E.S.D. BLOCK'/) 9999 IT=IT+1 RETURN END SUBROUTINE ZMENA(L,X,SX,K,IT,AT,IP) COMMON/UNITS/ W,M40,M50 COMMON/REFL/ YO,YC,A,B,DY,WDYQ,NULOVA,SINTHL,RHO,WT,DER(9708), 1 DC(9708),DYP,IQ COMMON NPS(9),NAM(9),NPSS(10),NAMM,NP,NZ,TLUM,Q,KI(9708),PS(9708), 1 AM(75000) EQUIVALENCE (NPSS(10),NPSM) LOGICAL NULOVA INTEGER W CHARACTER*8 AT,ATM CHARACTER*12 L,LM GO TO(1000,2000,6000),K 1000 IP=0 RXS=0. RXM=0. GO TO 9999 2000 IP=IP+1 CALL PIS(IT) DX=PS(IP) IF(KI(IP).EQ.0.AND.DX.EQ.0.) GO TO 5000 SX=DER(IP) RX=ABS(DX/SX) IF(RX.LE.RXM) GO TO 4000 RXM=RX ATM=AT LM=L 4000 RXS=RXS+RX XNEW=X+TLUM*DX WRITE(W,100) L,X,DX,XNEW,SX,RX,KI(IP) X=XNEW GO TO 9999 5000 WRITE(W,101) L,X,X GO TO 9999 6000 CALL NEWPG RXS=RXS/FLOAT(NPSM) WRITE(W,102) RXS,RXM,ATM,LM 100 FORMAT(1X,A12,1X,4F13.6,F10.2,I7) 101 FORMAT(1X,A12,1X, F13.6,13X,F13.6) 102 FORMAT('0AVERAGE CHANGE/E.S.D. IS',F7.2/'0MAXIMAL CHANGE/E.S.D. IS 1',F7.2,' FOR ',A8,1X,A12) 9999 RETURN END SUBROUTINE SMI(A,N,ISING,IBLOCK) DIMENSION A(1),H(200),D(200) COMMON/UNITS/ W,M40,M50 COMMON/COVC/ ICOV INTEGER W M=(N*(N+1))/2 ISING=0 IF(N.GT.1) GO TO 100 A(1)=1./A(1) RETURN 100 IJ=0 DO 500I=1,N IJ=IJ+I P=A(IJ) IF(P.LE.0.) GO TO 1010 500 D(I)=1./SQRT(P) CALL ZNORM(A,D,N) K=N+1 1000 K=K-1 IF(K.LE.0) GO TO 9000 P=A(1) IF(P.GT.0.)GO TO 1100 1010 ISING=1 RETURN 1100 P=1./P II=1 DO 3000I=2,N M=II II=II+I Q=A(M+1) HI=Q*P IF(I.LE.K) HI=-HI H(I)=HI M2=M+2 DO 2000IJ=M2,II IJI=IJ-I IJM=IJ-M 2000 A(IJI)=A(IJ)+Q*H(IJM) 3000 CONTINUE M=M-1 A(II)=P DO 4000I=2,N IJM=M+I 4000 A(IJM)=H(I) GO TO 1000 9000 CALL ZNORM(A,D,N) IF(ICOV.EQ.0) GO TO 9999 IJ=0 DO 9001I=1,N IJ=IJ+I P=A(IJ) 9001 D(I)=1./SQRT(P) CALL ZNORM(A,D,N) IC=0 K=0 DO 9020J=1,N DO 9010I=1,J K=K+1 IF(I.EQ.J) GO TO 9010 IF(ABS(A(K)).LE..9) GO TO 9010 IC=IC+1 IF(MOD(IC,50).NE.1) GO TO 9005 CALL NEWPG WRITE(W,101) IBLOCK 101 FORMAT('0ELEMENTS OF COVARIATION MATRIX (BLOCK #',I1, 1 ') LARGER THEN 0.9') 9005 CALL COVA(I,J,A(K),IBLOCK) 9010 CONTINUE 9020 CONTINUE DO 9030I=1,N 9030 D(I)=1./D(I) CALL ZNORM(A,D,N) 9999 RETURN END SUBROUTINE ZNORM(A,D,N) DIMENSION A(1),D(1) K=0 DO 2000J=1,N DJ=D(J) DO 1000I=1,J K=K+1 1000 A(K)=A(K)*DJ*D(I) 2000 CONTINUE RETURN END SUBROUTINE NASOB(A,P,X,N) DIMENSION A(1),P(1),X(1) IJP=0 DO 2000I=1,N XI=0.0 IJ=IJP ID=1 DO 1000J=1,N IJ=IJ+ID XI=XI+A(IJ)*P(J) IF(J.GE.I) ID=J 1000 CONTINUE X(I)=XI 2000 IJP=IJP+I RETURN END SUBROUTINE ZMDER(DER,RM,N) DIMENSION DER(1),RM(1),DR(6) K=0 DO 2000I=1,N DR(I)=0. DO 1000J=1,N K=K+1 1000 DR(I)=DR(I)+RM(K)*DER(J) 2000 CONTINUE DO 3000I=1,N 3000 DER(I)=DR(I) RETURN END SUBROUTINE ZMZPET(PS,DER,RM,N) DIMENSION PS(1),DER(1),RM(1),P(6),DR(6) DO 2000J=1,N K=J P(J)=0. DR(J)=0. DO 1000I=1,N P(J)=P(J)+RM(K)*PS(I) DR(J)=DR(J)+(RM(K)*DER(I))**2 1000 K=K+N 2000 CONTINUE DO 3000I=1,N PS(I)=P(I) 3000 DER(I)=SQRT(DR(I)) RETURN END FUNCTION BESJ(X,N,D,K) M=IABS(N) Z=X*.5 FN=1. IF(M.LE.1) GO TO 2000 DO 1000I=2,M 1000 FN=FN*FLOAT(I) 2000 A=Z**M/FN Z=Z**2 S=A DO 3000I=1,100 A=A*Z/FLOAT(I*(M+I)) IF(K.NE.1) A=-A S=S+A IF(ABS(A).LT.D) GO TO 4000 3000 CONTINUE 4000 IF(MOD(M,2).NE.0.AND.N.LT.0.AND.K.NE.1) S=-S BESJ=S RETURN END SUBROUTINE SMI3(A,AI,D) DIMENSION A(3,3),AI(3,3) DO 1I=1,3 K=1+MOD(I+1,3) M=1+MOD(I+3,3) DO 1J=1,3 L=1+MOD(J+1,3) N=1+MOD(J+3,3) 1 AI(I,J)=A(L,K)*A(N,M)-A(L,M)*A(N,K) D=0.0 DO 2I=1,3 2 D=D+A(I,1)*AI(1,I) DO 3I=1,3 DO 3J=1,3 3 AI(I,J)=AI(I,J)/D RETURN END SUBROUTINE SETAMP COMMON/KLICE/ NA,NS,NF,CSYM,VYH,NC,NOWR,IWQ,YOMIN,YOMAX,BLKOEF, 1 CENTR,NCYKL,SNLMN,SNLMX,KIM,DIF,KIW,LITE,IFSQ COMMON/CM4/ ATOM(100),AI(100),ISF(100),ITF(100),X(3,100), 1 BETA(6,100),SC(6),KMODX(100),UX(3,4,100),UY(3,4,100) 2 ,PHF(100),QCNTM(100),KMODS(100),A0(100),AX(4,100), 3 AY(4,100),KMODB(100),BX(6,4,100),BY(6,4,100), 4 KFX(100),KFB(100) COMMON/METRICS/ G(9),GI(9),TR(9),VOL COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10), 1 UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10), 2 PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10), 3 SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100), 4 NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10), 5 UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10) CHARACTER*8 ATOM,MOLNAME DIMENSION A(3) DO 2000I=1,NMOL K=KMODXM(I) IF(K.LE.0) GO TO 2000 DO 1000J=1,K A(1)=0. A(2)= URX(3,J,I)*VOL A(3)=-URX(2,J,I)*VOL CALL MULTM(GI,A,DURDX(1,J,I),3,3,1) A(2)= URY(3,J,I)*VOL A(3)=-URY(2,J,I)*VOL CALL MULTM(GI,A,DURDX(4,J,I),3,3,1) A(1)=-URX(3,J,I)*VOL A(2)=0. A(3)= URX(1,J,I)*VOL CALL MULTM(GI,A,DURDX(7,J,I),3,3,1) A(1)=-URY(3,J,I)*VOL A(3)= URY(1,J,I)*VOL CALL MULTM(GI,A,DURDX(10,J,I),3,3,1) A(1)= URX(2,J,I)*VOL A(2)=-URX(1,J,I)*VOL A(3)=0. CALL MULTM(GI,A,DURDX(13,J,I),3,3,1) A(1)= URY(2,J,I)*VOL A(2)=-URY(1,J,I)*VOL 1000 CALL MULTM(GI,A,DURDX(16,J,I),3,3,1) 2000 CONTINUE DO 9000I=1,NA IF(KMODX(I).GE.0) GO TO 9000 N=-KMODX(I) K=KMODXM(N) IF(K.LE.0) GO TO 9000 PHF(I)=PHFM(N) X1=X(1,I)-XM(1,N) X2=X(2,I)-XM(2,N) X3=X(3,I)-XM(3,N) DO 8000J=1,K DO 3000L=1,3 UX(L,J,I)=UTX(L,J,N) 3000 UY(L,J,I)=UTY(L,J,N) A(1)=-X2*URX(3,J,N)+X3*URX(2,J,N) A(2)=-X3*URX(1,J,N)+X1*URX(3,J,N) A(3)=-X1*URX(2,J,N)+X2*URX(1,J,N) DO 4000L=1,3 4000 A(L)=A(L)*VOL CALL CULTM(GI,A,UX(1,J,I),3,3,1) A(1)=-X2*URY(3,J,N)+X3*URY(2,J,N) A(2)=-X3*URY(1,J,N)+X1*URY(3,J,N) A(3)=-X1*URY(2,J,N)+X2*URY(1,J,N) DO 5000L=1,3 5000 A(L)=A(L)*VOL 8000 CALL CULTM(GI,A,UY(1,J,I),3,3,1) A(1)=0. A(2)=-X3*VOL A(3)= X2*VOL CALL MULTM(GI,A,DURDR(1,I),3,3,1) A(1)= X3*VOL A(2)=0. A(3)=-X1*VOL CALL MULTM(GI,A,DURDR(4,I),3,3,1) A(1)=-X2*VOL A(2)= X1*VOL A(3)=0. CALL MULTM(GI,A,DURDR(7,I),3,3,1) 9000 CONTINUE RETURN END SUBROUTINE CULTM(A,B,C,N1,N2,N3) DIMENSION A(1),B(1),C(1) DO 3000I=1,N1 JKP=1 IK=I DO 2000K=1,N3 IJ=I JK=JKP P=C(IK) DO 1000J=1,N2 P=P+A(IJ)*B(JK) IJ=IJ+N1 1000 JK=JK+1 C(IK)=P IK=IK+N1 2000 JKP=JKP+N2 3000 CONTINUE RETURN END SUBROUTINE MULTM(A,B,C,N1,N2,N3) DIMENSION A(1),B(1),C(1) DO 3000I=1,N1 JKP=1 IK=I DO 2000K=1,N3 IJ=I JK=JKP P=0. DO 1000J=1,N2 P=P+A(IJ)*B(JK) IJ=IJ+N1 1000 JK=JK+1 C(IK)=P IK=IK+N1 2000 JKP=JKP+N2 3000 CONTINUE RETURN END SUBROUTINE INDEX(J,K,L) K=J L=J IF(J.LE.3) RETURN K=J/3 L=(J+4)/3 RETURN END SUBROUTINE COVA(IA,IB,COV,IBLOCK) COMMON/KLICE/ NA,NS,NF,CSYM,VYH,NC,NOWR,IWQ,YOMIN,YOMAX,BLKOEF, 1 CENTR,NCYKL,SNLMN,SNLMX,KIM,DIF,KIW,LITE,IFSQ COMMON/UNITS/ W,M40,M50 COMMON/CM4/ ATOM(100),AI(100),ISF(100),ITF(100),X(3,100), 1 BETA(6,100),SC(6),KMODX(100),UX(3,4,100),UY(3,4,100) 2 ,PHF(100),QCNTM(100),KMODS(100),A0(100),AX(4,100), 3 AY(4,100),KMODB(100),BX(6,4,100),BY(6,4,100), 4 KFX(100),KFB(100) COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10), 1 UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10), 2 PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10), 3 SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100), 4 NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10), 5 UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10) COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10), 1 BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ, 2 EC(12) COMMON/LABEL/ L1,L2(3),L5,L7,L8,L9,L10 COMMON NPS(9),NAM(9),NPSS(10),NAMM,NP,NZ,TLUM,Q,KI(6883),PS(6883), 1 AM(75000) INTEGER W LOGICAL CSYM CHARACTER*1 CU,CV,CW,CP CHARACTER*4 CQ,S4 CHARACTER*8 MOLNAME,ATOM,S8,ATO CHARACTER*12 L1,L2,L5,L7,L8,L9,L10,LI CHARACTER*85 IVEN DIMENSION IP(2) DATA CU,CV,CW/'U','V','W'/,S8/' '/,S4/' '/ ENCODE(85,100,IVEN) IA,IB,COV 100 FORMAT(1X,'(',I3,',',I3,')',F7.3,' COVARIATION BETWEEN',48X) J=38 K=0 DO 1000I=1,NPMK IF(KI(I).NE.IBLOCK) GO TO 1000 K=K+1 IF(K.EQ.IA) IP(1)=I IF(K-IB) 1000,500,1500 500 IP(2)=I 1000 CONTINUE 1500 DO 3000I=1,2 ATO=S8 CQ=S4 K=IP(I) IF(K.GT.6) GO TO 2100 ENCODE(12,102,LI) K GO TO 2900 2100 IF(K.GT.18) GO TO 2200 K=K-6 IF(K.NE.1.AND.K.NE.7) GO TO 2140 IF(IEXT.NE.1) GO TO 2140 IF(K.EQ.7) GO TO 2120 LI=L7 GO TO 2900 2120 LI=L8 GO TO 2900 2140 IF(K.GT.6) GO TO 2150 CALL INDEX(K,II,JJ) ENCODE(12,110,LI) II,JJ GO TO 2900 2150 CALL INDEX(K-6,II,JJ) ENCODE(12,111,LI) II,JJ GO TO 2900 2200 IF(K.GE.NPMP) GO TO 2700 II=(K-19)/92+1 ATO=ATOM(II) K=MOD(K-19,92)+1 IF(K.GT.1) GO TO 2300 LI=L1 GO TO 2900 2300 IF(K.GT.4) GO TO 2400 LI=L2(K-1) GO TO 2900 2400 IF(K.GT.10) GO TO 2500 CALL INDEX(K-4,II,JJ) IF(LITE.EQ.1) GO TO 2450 ENCODE(12,112,LI) II,JJ GO TO 2900 2450 ENCODE(12,117,LI) II,JJ GO TO 2900 2500 IF(K.GT.19) GO TO 2530 IF(K.GT.11) GO TO 2510 LI=L10 GO TO 2900 2510 II=(K-12)/2+1 JJ=MOD(K-12,2)+1 IF(JJ.GT.1) GO TO 2520 ENCODE(12,115,LI) II GO TO 2900 2520 ENCODE(12,116,LI) II GO TO 2900 2530 IF(K.GT.43) GO TO 2590 II=(K-20)/6+1 JJ=MOD(K-20,6)+1 IF(JJ.GT.3) GO TO 2580 ENCODE(12,113,LI) CU,JJ,II GO TO 2900 2580 ENCODE(12,114,LI) CU,JJ-3,II GO TO 2900 2590 IF(K.GT.91) GO TO 2600 II=(K-44)/12+1 JJ=MOD(K-44,12)+1 IF(JJ.LE.6) THEN ENCODE(12,113,LI) 'T',JJ,II ELSE ENCODE(12,113,LI) 'T',JJ-6,II ENDIF GO TO 2900 2600 LI=L9 GO TO 2900 2700 II=(K-NPMP)/49+1 ATO=MOLNAME(II) K=MOD(K-NPMP,49)+1 IF(K.GT.48) GO TO 2800 CP=CV IF(K.LE.24) GO TO 2720 K=K-24 CP=CW 2720 II=(K-1)/6+1 JJ=MOD(K-1,6)+1 IF(JJ.GT.3) GO TO 2730 ENCODE(12,113,LI) CP,JJ,II GO TO 2900 2730 ENCODE(12,114,LI) CP,JJ-3,II GO TO 2900 2800 LI=L9 2900 IF(I.EQ.1) GO TO 2910 IVEN=IVEN(1:J)//' AND ' J=J+5 2910 K=IDEL(LI) IF(K.NE.0) IVEN=IVEN(1:J)//LI(1:K) J=J+K+1 K=IDEL(ATO) IF(K.NE.0) IVEN=IVEN(1:J)//ATO(1:K) J=J+K+1 J=IDEL(IVEN) 3000 CONTINUE WRITE(W,101) IVEN RETURN 101 FORMAT(A85) 102 FORMAT('SCALE-',I1) 110 FORMAT('RHO(',I1,',',I1,') ') 111 FORMAT('G(',I1,',',I1,') ') 112 FORMAT('BETA(',I1,',',I1,') ') 113 FORMAT(A1,I1,'*SIN(',I1,'*X4)') 114 FORMAT(A1,I1,'*COS(',I1,'*X4)') 115 FORMAT('A*SIN(',I1,'*X4) ') 116 FORMAT('A*COS(',I1,'*X4) ') 117 FORMAT('U(',I1,',',I1,') ') END FUNCTION IDEL(L) CHARACTER*(*) L IL=LEN(L) 500 IF(L(IL:IL).NE.' ') GO TO 1000 IL=IL-1 IF(IL.GT.0) GO TO 500 1000 IDEL=IL RETURN END SUBROUTINE CENTER(I) COMMON/KLICE/ NA,NS,NF,CSYM,VYH,NC,NOWR,IWQ,YOMIN,YOMAX,BLKOEF, 1 CENTR,NCYKL,SNLMN,SNLMX,KIM,DIF,KIW,LITE,IFSQ COMMON/ATOMY/ FF(32,15),FX(15),FY(15),AW(15) COMMON/CM4/ ATOM(100),AI(100),ISF(100),ITF(100),X(3,100), 1 BETA(6,100),SC(6),KMODX(100),UX(3,4,100),UY(3,4,100) 2 ,PHF(100),QCNTM(100),KMODS(100),A0(100),AX(4,100), 3 AY(4,100),KMODB(100),BX(6,4,100),BY(6,4,100), 4 KFX(100),KFB(100) COMMON/SYMM/ IS(2,4,24),TS(3,24),TAU(24),FSM(24) COMMON/METRICS/ G(9),GI(9),TR(9),VOL COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10), 1 UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10), 2 PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10), 3 SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100), 4 NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10), 5 UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10) COMMON/UNITS/ W,M40,M50 CHARACTER*8 ATOM,MOLNAME INTEGER W LOGICAL CSYM DIMENSION XP(3),A(3),TI(6),CP(3),TRM(9) DO 910J=1,6 910 TI(J)=0. IF(KIW.LE.0) GO TO 3500 FN=0. DO 900J=1,3 900 CP(J)=0. DO 3000J=1,NPG(I) DO 2000K=1,NA IF(KMODX(K).NE.-I) GO TO 2000 AIW=AI(K)*AW(ISF(K)) IF(J.GT.1) GO TO 1000 DO 950L=1,3 950 XP(L)=X(L,K) GO TO 1100 1000 M=KPG(J,I) MM=IABS(M) DO 1010L=1,3 XP(L)=0. 1010 A(L)=X(L,K)-XM(L,I) DO 1013L=1,3 DO 1012LL=1,2 IND=IS(LL,L,MM) IF(IND.EQ.0) GO TO 1012 Z=SIGN(1.,FLOAT(IND)) IND=IABS(IND) XP(L)=XP(L)+Z*A(IND) 1012 CONTINUE 1013 CONTINUE IF(M.GE.0) GO TO 1015 DO 1014L=1,3 1014 XP(L)=-XP(L) 1015 DO 1020L=1,3 1020 XP(L)=XP(L)+XM(L,I) 1100 FN=FN+AIW DO 1500L=1,3 1500 CP(L)=CP(L)+AIW*XP(L) 2000 CONTINUE 3000 CONTINUE DO 3100J=1,3 3100 XM(J,I)=CP(J)/FN 3500 DO 6000J=1,NPG(I) DO 5000K=1,NA IF(KMODX(K).NE.-I) GO TO 5000 AIW=AI(K)*AW(ISF(K)) IF(J.GT.1) GO TO 4000 DO 3550L=1,3 3550 XP(L)=X(L,K)-XM(L,I) GO TO 4100 4000 M=KPG(J,I) MM=IABS(M) DO 4010L=1,3 XP(L)=0. 4010 A(L)=X(L,K)-XM(L,I) DO 4013L=1,3 DO 4012LL=1,2 IND=IS(LL,L,MM) IF(IND.EQ.0) GO TO 4012 Z=SIGN(1.,FLOAT(IND)) IND=IABS(IND) XP(L)=XP(L)+Z*A(IND) 4012 CONTINUE 4013 CONTINUE IF(M.GE.0) GO TO 4100 DO 4014L=1,3 4014 XP(L)=-XP(L) 4100 CALL MULTM(TR,XP,A,3,3,1) TI(1)=TI(1)+AIW*(A(2)**2+A(3)**2) TI(2)=TI(2)-AIW*A(1)*A(2) TI(3)=TI(3)+AIW*(A(1)**2+A(3)**2) TI(4)=TI(4)-AIW*A(1)*A(3) TI(5)=TI(5)-AIW*A(2)*A(3) TI(6)=TI(6)+AIW*(A(1)**2+A(2)**2) 5000 CONTINUE 6000 CONTINUE CALL QL(TI,TRM,A,3) CALL TRMAT(TRM,TRMA(1,I),3) RETURN END SUBROUTINE INTER COMMON/UNITS/ W,M40,M50 COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10), 1 BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ, 2 EC(12) COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10), 1 UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10), 2 PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10), 3 SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100), 4 NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10), 5 UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10) COMMON/METRICS/ G(9),GI(9),TR(9),VOL DIMENSION X(3),Y(3),Z(3),SV(3),SM(9),S1(3),S2(3),TRMAP(9),SX(3), 1 SY(3),SZ(3) INTEGER W CHARACTER*8 MOLNAME IF(NMOL.LE.0) GO TO 9999 CALL NEWPG WRITE(W,100) 100 FORMAT('0INTERPRETATION OF DISPLACEMENT VECTORS'/39('=')//) WRITE(W,106) 106 FORMAT('0COMPONENTS OF THE AMPLITUDE OF THE DISPLACEMENT VECTORS A 1RE EXPRESSED AS :'/'0X(I) .... I=1,2,3 THE FRACTIONAL COORDINATE 2 SYSTEM (NORMALISED TO BE THE MAXIMUM COORDINATE EQUAL 1) '/'0Y(I) 3 .... I=1,2,3 THE CARTESIAN COORDINATE CRYSTAL SYSTEM (NORMALISE 4 TO BE LENGTH EQUAL 1)'/'0Z(I) .... I=1,2,3 THE CARTESIAN COORDI 5NATE TENSOR INERTIA SYSTEM (NORMALISE TO BE LENGTH EQUAL 1)') WRITE(W,107) 107 FORMAT(/' MAGNITUDE OF DISPLACEMENT VECTORS IS DONE IN ANGSTROMS F 1OR TRANSLATION AND IN DEGREES FOR ROTATIONAL DISPLACEMENT'/) DO 9000I=1,NMOL WRITE(W,101) MOLNAME(I) 101 FORMAT(//' MOLECULE : ',A8/20('=')/) WRITE(W,102) 102 FORMAT(/' TRANSLATION DISPLACEMENT VECTORS'/) WRITE(W,108) 108 FORMAT(' WAVE MAGNITUDE',4X,'X(1)',4X,'X(2)',4X,'X(3)', 14X,'Y(1)',4X,'Y(2)',4X,'Y(3)',4X,'Z(1)',4X,'Z(2)',4X,'Z(3)'/ 2 1X,91('-')) CALL MULTM(TRMA(1,I),TR,TRMAP,3,3,3) DO 1300K=1,KMODXM(I) CALL NORM(UTX(1,K,I),X,DX) IF(DX.EQ.0.) GO TO 1100 DO 500L=1,3 500 SX(L)=SUTX(L,K,I)/DX CALL MULTM(TR,UTX(1,K,I),Y,3,3,1) D=SQRT(Y(1)**2+Y(2)**2+Y(3)**2) DO 1000L=1,3 1000 Y(L)=Y(L)/D CALL DODEJ(TR,D,Y,SV,SM) IP=1 DO 1001L=1,3 CALL NASOB(UTXM(1,K,I),SM(IP),S1,3) CALL MULTM(SM(IP),S1,SY(L),1,3,1) 1001 IP=IP+3 CALL NASOB(UTXM(1,K,I),SV,S1,3) CALL MULTM(SV,S1,SD,1,3,1) CALL MULTM(TRMAP,UTX(1,K,I),Z,3,3,1) DO 1010L=1,3 1010 Z(L)=Z(L)/D CALL DODEJ(TRMAP,D,Z,SV,SM) IP=1 DO 1011L=1,3 CALL NASOB(UTXM(1,K,I),SM(IP),S1,3) CALL MULTM(SM(IP),S1,SZ(L),1,3,1) 1011 IP=IP+3 DO 1020L=1,3 SY(L)=SQRT(SY(L)) 1020 SZ(L)=SQRT(SZ(L)) SD=SQRT(SD) WRITE(W,103) K,D,X,Y,Z,SD,SX,SY,SZ 103 FORMAT(' SIN(',I1,'X4)',F10.5,9F8.3/9X,F10.5,9F8.3) 1100 CALL NORM(UTY(1,K,I),X,DX) IF(DX.EQ.0.) GO TO 1300 DO 1110L=1,3 1110 SX(L)=SUTY(L,K,I)/DX CALL MULTM(TR,UTY(1,K,I),Y,3,3,1) D=SQRT(Y(1)**2+Y(2)**2+Y(3)**2) DO 1120L=1,3 1120 Y(L)=Y(L)/D CALL DODEJ(TR,D,Y,SV,SM) IP=1 DO 1021L=1,3 CALL NASOB(UTYM(1,K,I),SM(IP),S1,3) CALL MULTM(SM(IP),S1,SY(L),1,3,1) 1021 IP=IP+3 CALL NASOB(UTYM(1,K,I),SV,S1,3) CALL MULTM(SV,S1,SD,1,3,1) CALL MULTM(TRMAP,UTY(1,K,I),Z,3,3,1) DO 1130L=1,3 1130 Z(L)=Z(L)/D CALL DODEJ(TRMAP,D,Z,SV,SM) IP=1 DO 1131L=1,3 CALL NASOB(UTYM(1,K,I),SM(IP),S1,3) CALL MULTM(SM(IP),S1,SZ(L),1,3,1) 1131 IP=IP+3 DO 1140L=1,3 SY(L)=SQRT(SY(L)) 1140 SZ(L)=SQRT(SZ(L)) SD=SQRT(SD) WRITE(W,110) K,D,X,Y,Z,SD,SX,SY,SZ 110 FORMAT(' COS(',I1,'X4)',F10.5,9F8.3/9X,F10.5,9F8.3) 1300 CONTINUE WRITE(W,104) 104 FORMAT(/' ROTATIONAL DISPLACEMENT VECTORS'/) WRITE(W,108) DO 2300K=1,KMODXM(I) CALL NORM(URX(1,K,I),X,DX) IF(DX.EQ.0.) GO TO 2100 DO 1500L=1,3 1500 SX(L)=SURX(L,K,I)/DX CALL MULTM(TR,URX(1,K,I),Y,3,3,1) D=SQRT(Y(1)**2+Y(2)**2+Y(3)**2) DO 2000L=1,3 2000 Y(L)=Y(L)/D CALL DODEJ(TR,D,Y,SV,SM) IP=1 DO 2001L=1,3 CALL NASOB(URXM(1,K,I),SM(IP),S1,3) CALL MULTM(SM(IP),S1,SY(L),1,3,1) 2001 IP=IP+3 CALL NASOB(URXM(1,K,I),SV,S1,3) CALL MULTM(SV,S1,SD,1,3,1) CALL MULTM(TRMAP,URX(1,K,I),Z,3,3,1) DO 2010L=1,3 2010 Z(L)=Z(L)/D CALL DODEJ(TRMAP,D,Z,SV,SM) IP=1 DO 2011L=1,3 CALL NASOB(URXM(1,K,I),SM(IP),S1,3) CALL MULTM(SM(IP),S1,SZ(L),1,3,1) 2011 IP=IP+3 DO 2020L=1,3 SY(L)=SQRT(SY(L)) 2020 SZ(L)=SQRT(SZ(L)) SD=SQRT(SD)/TORAD D=D/TORAD WRITE(W,109) K,D,X,Y,Z,SD,SX,SY,SZ 109 FORMAT(' SIN(',I1,'X4)',F10.2,9F8.3/9X,F10.2,9F8.3) 2100 CALL NORM(URY(1,K,I),X,DX) IF(DX.EQ.0.) GO TO 2300 DO 2110L=1,3 2110 SX(L)=SURY(L,K,I)/DX CALL MULTM(TR,URY(1,K,I),Y,3,3,1) D=SQRT(Y(1)**2+Y(2)**2+Y(3)**2) DO 2120L=1,3 2120 Y(L)=Y(L)/D CALL DODEJ(TR,D,Y,SV,SM) IP=1 DO 2121L=1,3 CALL NASOB(URYM(1,K,I),SM(IP),S1,3) CALL MULTM(SM(IP),S1,SY(L),1,3,1) 2121 IP=IP+3 CALL NASOB(URYM(1,K,I),SV,S1,3) CALL MULTM(SV,S1,SD,1,3,1) CALL MULTM(TRMAP,URY(1,K,I),Z,3,3,1) DO 2130L=1,3 2130 Z(L)=Z(L)/D CALL DODEJ(TRMAP,D,Z,SV,SM) IP=1 DO 2131L=1,3 CALL NASOB(URYM(1,K,I),SM(IP),S1,3) CALL MULTM(SM(IP),S1,SZ(L),1,3,1) 2131 IP=IP+3 DO 2140L=1,3 SY(L)=SQRT(SY(L)) 2140 SZ(L)=SQRT(SZ(L)) SD=SQRT(SD)/TORAD D=D/TORAD WRITE(W,105) K,D,X,Y,Z,SD,SX,SY,SZ 105 FORMAT(' COS(',I1,'X4)',F10.2,9F8.3/9X,F10.2,9F8.3) 2300 CONTINUE 9000 CONTINUE 9999 RETURN END SUBROUTINE NORM(X,Y,D) DIMENSION X(3),Y(3) D=ABS(X(1)) DO 1000 I=2,3 IF(D.LT.ABS(X(I))) D=ABS(X(I)) 1000 CONTINUE IF(D.EQ.0.) RETURN DO 2000 I=1,3 2000 Y(I)=X(I)/D RETURN END SUBROUTINE DODEJ(T,D,Z,SV,SM) DIMENSION T(9),Z(3),SV(3),SM(9),SMP(9) DO 1000I=1,9 1000 SMP(I)=T(I) CALL MULTM(Z,T,SV,1,3,3) DO 1500I=1,3 1500 SV(I)=-SV(I) CALL CULTM(SV,Z,SMP,3,1,3) DO 2000I=1,9 2000 SMP(I)=SMP(I)/D CALL TRMAT(SMP,SM,3) RETURN END SUBROUTINE SKOVEJ(AM,Q,IBLOCK,KI) COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10), 1 UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10), 2 PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10), 3 SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100), 4 NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10), 5 UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10) DIMENSION AM(1),KI(1) CHARACTER*8 MOLNAME NK=0 DO 500I=1,NPMP-1 IF(KI(I).EQ.IBLOCK) NK=NK+1 500 CONTINUE IP=NPMP DO 5000I=1,NMOL DO 3000K=1,4 IF(KI(IP ).NE.IBLOCK.AND.KI(IP+1).NE.IBLOCK 1 .AND.KI(IP+2).NE.IBLOCK) GO TO 1000 CALL KTERE(KI,IBLOCK,NK,IP,I1,I2,I3) CALL ZAPLN(AM,I1,I2,I3,Q,UTXM(1,K,I)) 1000 IF(KI(IP+3).NE.IBLOCK.AND.KI(IP+4).NE.IBLOCK 1 .AND.KI(IP+5).NE.IBLOCK) GO TO 2000 CALL KTERE(KI,IBLOCK,NK,IP+3,I1,I2,I3) CALL ZAPLN(AM,I1,I2,I3,Q,UTYM(1,K,I)) 2000 IF(KI(IP+24).NE.IBLOCK.AND.KI(IP+25).NE.IBLOCK 1 .AND.KI(IP+26).NE.IBLOCK) GO TO 2500 CALL KTERE(KI,IBLOCK,NK,IP+24,I1,I2,I3) CALL ZAPLN(AM,I1,I2,I3,Q,URXM(1,K,I)) 2500 IF(KI(IP+27).NE.IBLOCK.AND.KI(IP+28).NE.IBLOCK 1 .AND.KI(IP+29).NE.IBLOCK) GO TO 3000 CALL KTERE(KI,IBLOCK,NK,IP+27,I1,I2,I3) CALL ZAPLN(AM,I1,I2,I3,Q,URYM(1,K,I)) 3000 IP=IP+6 5000 IP=IP+25 RETURN END SUBROUTINE KTERE(KI,IBLOCK,NK,IP,I1,I2,I3) DIMENSION KI(1) I1=0 IF(KI(IP).NE.IBLOCK) GO TO 1000 NK=NK+1 I1=NK 1000 I2=0 IF(KI(IP+1).NE.IBLOCK) GO TO 2000 NK=NK+1 I2=NK 2000 I3=0 IF(KI(IP+2).NE.IBLOCK) GO TO 3000 NK=NK+1 I3=NK 3000 RETURN END SUBROUTINE ZAPLN(A,I1,I2,I3,Q,U) DIMENSION A(1),U(1) IMAT(II,JJ)=II+(JJ*(JJ-1))/2 IF(I1.EQ.0) GO TO 2000 U(1)=A(IMAT(I1,I1))*Q IF(I2.EQ.0) GO TO 1000 U(2)=A(IMAT(I1,I2))*Q 1000 IF(I3.EQ.0) GO TO 2000 U(4)=A(IMAT(I1,I3))*Q 2000 IF(I2.EQ.0) GO TO 3000 U(3)=A(IMAT(I2,I2))*Q IF(I3.EQ.0) GO TO 9999 U(5)=A(IMAT(I2,I3))*Q 3000 IF(I3.EQ.0) GO TO 9999 U(6)=A(IMAT(I3,I3))*Q 9999 RETURN END SUBROUTINE RECIP(CP,RCP) DIMENSION CP(6),RCP(6) AR=CP(1) BR=CP(2) CR=CP(3) COSAR=CP(4) COSBR=CP(5) COSGR=CP(6) VR=AR*BR*CR*SQRT(1.+2.*COSAR*COSBR*COSGR- 1 COSAR**2-COSBR**2-COSGR**2) SINAR=SQRT(1.-COSAR**2) SINBR=SQRT(1.-COSBR**2) SINGR=SQRT(1.-COSGR**2) A=BR*CR*SINAR/VR B=CR*AR*SINBR/VR C=AR*BR*SINGR/VR COSA=(COSBR*COSGR-COSAR)/(SINBR*SINGR) COSB=(COSGR*COSAR-COSBR)/(SINGR*SINAR) COSC=(COSAR*COSBR-COSGR)/(SINAR*SINBR) RCP(1)=A RCP(2)=B RCP(3)=C RCP(4)=COSA RCP(5)=COSB RCP(6)=COSC RETURN END SUBROUTINE RESETX(X) COMMON/UNITS/ W,M40,M50 INTEGER W DIMENSION X(3,1) DATA IFIRST/0/ IF(IFIRST.EQ.0) WRITE(W,100) 100 FORMAT('0YOU USED A DUMMY SUBROUTINE RESETX') IFIRST=1 RETURN END SUBROUTINE RESETA(A0,AX,AY) COMMON/UNITS/ W,M40,M50 INTEGER W DIMENSION A0(1),AX(4,1),AY(4,1) DATA IFIRST/0/ IF(IFIRST.EQ.0) WRITE(W,100) 100 FORMAT('0YOU USED A DUMMY SUBROUTINE RESETA') IFIRST=1 RETURN SUBROUTINE RESETC(BX,BY) COMMON/UNITS/ W,M40,M50 INTEGER W DIMENSION BX(6,6,1),BY(6,4,1) DATA IFIRST/0/ IF(IFIRST.EQ.0) WRITE(W,100) 100 FORMAT('0YOU USED A DUMMY SUBROUTINE RESETC') IFIRST=1 RETURN END SUBROUTINE RESETB(BETA) COMMON/UNITS/ W,M40,M50 INTEGER W DIMENSION BETA(6,1) DATA IFIRST/0/ IF(IFIRST.EQ.0) WRITE(W,100) 100 FORMAT('0YOU USED A DUMMY SUBROUTINE RESETB') IFIRST=1 RETURN END SUBROUTINE RESETD(UX,UY) COMMON/UNITS/ W,M40,M50 INTEGER W DIMENSION UX(3,4,1),UY(3,4,1) DATA IFIRST/0/ IF(IFIRST.EQ.0) WRITE(W,100) 100 FORMAT('0YOU USED A DUMMY SUBROUTINE RESETD') IFIRST=1 RETURN END SUBROUTINE RESETM(UTX,UTY,URX,URY) COMMON/UNITS/ W,M40,M50 INTEGER W DIMENSION UTX(3,4,1),UTY(3,4,1),URX(3,4,1),URY(3,4,1) DATA IFIRST/0/ IF(IFIRST.EQ.0) WRITE(W,100) 100 FORMAT('0YOU USED A DUMMY SUBROUTINE RESETM') IFIRST=1 RETURN END SUBROUTINE DSETA(DER,NA) DIMENSION DER(92,NA) C C SUBROUTINE FOR SETTING OF DERIVATIONS OF ATOMIC PARAMETERS C C FOR EACH ATOM WE HAVE 92 PARAMETERS : C C #1 AI - OCCUPATIONAL PARAMETER C #2-#4 X(I),I=1,3 - POSITIONAL PARAMETERS C #5-#10 BETA(I),I=1,6 - ANISOTROPIC TEMPERATURE PARAMETERS C IF ISOTROPIC ATOM ONLY #5 IS USED C AND #6-#10 ARE VOID C *********************************** C **** ATOMIC PARAMETERS FOR **** C **** OCCUPATIONAL MODULATION **** C *********************************** C #11 A0 - ABSOLUTE TERM C #12 AX(1) - SIN WAVE OF 1ST HARMONIC C #13 AY(1) - COS WAVE OF 1ST HARMINIC C #14 AX(2) - SIN WAVE OF 2ND HARMONIC C #15 AY(2) - COS WAVE OF 2ND HARMINIC C #16 AX(3) - SIN WAVE OF 3RD HARMONIC C #17 AY(3) - COS WAVE OF 3RD HARMINIC C #18 AX(4) - SIN WAVE OF 4TH HARMONIC C #19 AY(4) - COS WAVE OF 4TH HARMINIC C *********************************** C **** ATOMIC PARAMETERS FOR **** C **** POSITIONAL MODULATION **** C *********************************** C #20-#22 UX(I,1,1),I=1,3 - SIN WAVE OF 1ST HARMONIC C #23-#25 UY(I,1,1),I=1,3 - COS WAVE OF 1ST HARMINIC C #26-#28 UX(I,2,1),I=1,3 - SIN WAVE OF 2ND HARMONIC C #29-#31 UY(I,2,1),I=1,3 - COS WAVE OF 2ND HARMINIC C #32-#34 UX(I,3,1),I=1,3 - SIN WAVE OF 3RD HARMONIC C #35-#37 UY(I,3,1),I=1,3 - COS WAVE OF 3RD HARMINIC C #38-#40 UX(I,4,1),I=1,3 - SIN WAVE OF 4TH HARMONIC C #41-#43 UY(I,4,1),I=1,3 - COS WAVE OF 4TH HARMINIC C *********************************** C **** ATOMIC PARAMETERS FOR **** C **** TEMP.FACTOR MODULATION **** C *********************************** C #44-#49 BX(I,4,1),I=1,6 - SIN WAVE OF 1ST HARMONIC C #50-#55 BY(I,4,1),I=1,6 - COS WAVE OF 1ST HARMONIC C #56-#61 BX(I,4,2),I=1,6 - SIN WAVE OF 2ND HARMONIC C #62-#67 BY(I,4,2),I=1,6 - COS WAVE OF 2ND HARMONIC C #68-#73 BX(I,4,3),I=1,6 - SIN WAVE OF 3RD HARMONIC C #74-#79 BY(I,4,3),I=1,6 - COS WAVE OF 3RD HARMONIC C #80-#85 BX(I,4,4),I=1,6 - SIN WAVE OF 4TH HARMONIC C #86-#91 BY(I,4,4),I=1,6 - COS WAVE OF 4TH HARMONIC C #92 PH - PHASON FACTOR RETURN END SUBROUTINE DSETM(DER,NMOL) DIMENSION DER(49,NMOL) C C SUBROUTINE FOR SETTING OF DERIVATIONS OF MOLECULAR PARAMETERS C C FOR EACH MOLECULE WE HAVE 49 PARAMETERS : C C #1-#3 UTX(I,1,1),I=1,3 - TRANSLATIONAL SIN WAVE OF 1ST HARMONIC C #4-#6 UTY(I,1,1),I=1,3 - TRANSLATIONAL COS WAVE OF 1ST HARMONIC C #7-#10 UTX(I,2,1),I=1,3 - TRANSLATIONAL SIN WAVE OF 2ND HARMONIC C #11-#12 UTY(I,2,1),I=1,3 - TRANSLATIONAL COS WAVE OF 2ND HARMONIC C #13-#15 UTX(I,3,1),I=1,3 - TRANSLATIONAL SIN WAVE OF 3RD HARMONIC C #16-#18 UTY(I,3,1),I=1,3 - TRANSLATIONAL COS WAVE OF 3RD HARMONIC C #19-#21 UTX(I,4,1),I=1,3 - TRANSLATIONAL SIN WAVE OF 4TH HARMONIC C #25-#27 URX(I,1,1),I=1,3 - ROTATIONAL SIN WAVE OF 1ST HARMONIC C #28-#30 URY(I,1,1),I=1,3 - ROTATIONAL COS WAVE OF 1ST HARMONIC C #31-#33 URX(I,2,1),I=1,3 - ROTATIONAL SIN WAVE OF 2ND HARMONIC C #34-#36 URY(I,2,1),I=1,3 - ROTATIONAL COS WAVE OF 2ND HARMONIC C #37-#39 URX(I,3,1),I=1,3 - ROTATIONAL SIN WAVE OF 3RD HARMONIC C #40-#42 URY(I,3,1),I=1,3 - ROTATIONAL COS WAVE OF 3RD HARMONIC C #43-#45 URX(I,4,1),I=1,3 - ROTATIONAL SIN WAVE OF 4TH HARMONIC C #46-#48 URY(I,4,1),I=1,3 - ROTATIONAL COS WAVE OF 4TH HARMONIC C #49 PH - MOLECULAR PHASON FACTOR RETURN END SUBROUTINE TRMAT(A,B,N) C C .... SUBROUTINE FOR CALCULATION OF TRANSPOSE MATRIX C DIMENSION A(N,N),B(N,N) DO 2000J=1,N DO 1000I=1,N 1000 B(J,I)=A(I,J) 2000 CONTINUE RETURN END SUBROUTINE QL(A,Z,D,N) DIMENSION E(20),D(1),A(1),Z(1) DATA TOL,EPS/.0000001,.0000001/ IMD=1-2*MOD(N,2) NB=N*(N-1)+1 K=0 DO 20I=1,N DO 10J=1,N K=K+1 IF(I-J) 8,9,8 9 Z(K)=1. GO TO 10 8 Z(K)=0. 10 CONTINUE 20 CONTINUE IZ=N*(N+1)/2 L=N DO 900IP=3,N I=L L=L-1 H=0. IZ=IZ-I DO 100K=1,L IND=IZ+K F=A(IND) D(K)=F 100 H=H+F**2 IF(H.GT.TOL) GO TO 150 E(I)=0. H=0. GO TO 890 150 G=-SQRT(H) IF(F.LT.0.) G=-G E(I)=G H=H-F*G RH=1./H DL=F-G D(L)=DL IND=IZ+L A(IND)=DL F=0. IND=1 DO 300J=1,L G=0. JK=IND DO 200K=1,L G=G+A(JK )*D(K) ID=1 IF(J.LE.K) ID=K JK=JK +ID 200 CONTINUE G=G*RH E(J)=G F=F+G*D(J) IND=IND+J 300 CONTINUE HH=F*.5*RH JK=0 DO 500J=1,L F=D(J) G=E(J)-HH*F E(J)=G DO 400K=1,J JK=JK+1 400 A(JK)=A(JK)-F*E(K)-G*D(K) 500 CONTINUE 890 JK=IZ+I D(I)=A(JK) A(JK)=H 900 CONTINUE D(1)=A(1) E(2)=A(2) D(2)=A(3) B=0. F=0. DO 2000L=1,N LN=L*N-N IF(L.EQ.N) GO TO 1800 L1=L+1 H=EPS*(ABS(D(L))+ABS(E(L1))) IF(B.LT.H) B=H DO 1100M=L1,N IF(ABS(E(M)).LE.B) GO TO 1200 1100 CONTINUE M=N+1 1200 M=M-1 IF(M.EQ.L) GO TO 1800 ML=M-L MN=M*N 1250 G=D(L) EL=E(L1) DL=(D(L1)-G)/(2.*EL) H=ABS(DL)+SQRT(DL**2+1.) IF(DL.LT.0.) H=-H DL=EL/H D(L)=DL H=G-DL IF(L1.GT.N) GO TO 1350 DO 1300I=L1,N 1300 D(I)=D(I)-H 1350 F=F+H DL=D(M) C=1. S=0. IND=MN DO 1600JK=1,ML IND=IND-N I=M-JK I1=I+1 EL=E(I1) G=C*EL H=C*DL IF(ABS(DL).LT.ABS(EL)) GO TO 1400 C=EL/DL R=SQRT(C**2+1.) E(I+2)=S*DL*R S=C/R C=1./R GO TO 1500 1400 C=DL/EL R=SQRT(C**2+1.) E(I+2)=S*EL*R S=1./R C=C*S 1500 DL=C*D(I)-S*G D(I1 )=H+S*(C*G+S*D(I)) DO 1550K=1,N IZ=K+IND IZN=IZ-N H=Z(IZN) G=Z(IZ) Z(IZN) =C*H-S*G 1550 Z(IZ )=S*H+C*G 1600 CONTINUE EL=S*DL E(L1 )=EL D(L)=C*DL IF(ABS(EL).GT.B)GO TO 1250 1800 DL=D(L)+F JK=L-1 IF(JK.LT.1) GO TO 1855 IZ=LN DO 1850K=1,JK I=L-K IZ=IZ-N IF(DL.GE.D(I)) GO TO 1860 IMD=-IMD DO 1820J=1,N IND=IZ+J INZ=IND+N P=Z(IND) Z(IND)=Z(INZ) 1820 Z(INZ)=P 1850 D(I+1)=D(I) 1855 I=0 1860 D(I+1)=DL 2000 CONTINUE IF(IMD) 2010,2100,2100 2010 DO 2020I=1,N 2020 Z(I)=-Z(I) 2100 IZ=1 DO 2500I=3,N L=I-1 IZ=IZ+L IND=IZ+I H=A(IND) IF(H.EQ.0.) GO TO 2500 DO 2400J=1,NB,N JK=J-1 S=0. DO 2200K=1,L M=IZ+K IND=JK+K 2200 S=S+A(M)*Z(IND) S=S/H DO 2300K=1,L IND=JK+K M=IZ+K 2300 Z(IND)=Z(IND)-S*A(M) 2400 CONTINUE 2500 CONTINUE RETURN END SUBROUTINE BMOD(F,DF1,DF2,DF3,DB,D1,D2,M,ITYP) COMMON/STALY/ PI2 GO TO (1000,2000,3000),ITYP 1000 POM=EXP(DB) IF(M.EQ.0) THEN F=D1*POM+D2*(POM-1.)/DB-D1-D2+1. DF1=D1*POM+D2*(POM*DB-POM+1.)/DB**2 DF2=POM-1. DF3=(POM-1.)/DB-1. ELSE ARG=PI2*FLOAT(M)*.5 ALFA1=ARG*D1 CS1=COS(ALFA1) SN1=SIN(ALFA1) IF(D2.NE.0.) THEN ALFA2=ARG*D2 AQ=ALFA2**2 CS12=COS(ALFA1+ALFA2) SN12=SIN(ALFA1+ALFA2) BQ=DB**2 AB=ALFA2*DB FCIT=AB*POM*CS1-AB*CS12+BQ*POM*SN1-BQ*SN12 FJMEN=ALFA2*(AQ+BQ) FJQ=FJMEN**2 Z=FCIT/FJMEN F=D2*Z DF1=D2*((POM*CS1*(AB+ALFA2)-ALFA2*CS12+2.*DB*(POM*SN1-SN12))*FJMEN 1 -2.*AB*FCIT)/FJQ DF2=(-AB*ALFA2*POM*SN1+AB*ALFA2*SN12+BQ*POM*ALFA2*CS1 1 -BQ*ALFA2*CS12)/FJMEN DF3=Z+((AB*POM*CS1-AB*CS12+ALFA2*AB*SN12-BQ*ALFA2*CS12)*FJMEN 1 -(3.*AQ+BQ)*ALFA2*FCIT)/FJQ ELSE FCIT=SN1/ARG F=(POM-1.)*FCIT DF1=POM*FCIT DF2=(POM-1.)*CS1 DF3=0. ENDIF ENDIF RETURN 2000 IF(DB.NE.0.) THEN D=SQRT(ABS(DB)) POM=SIN(D)/D POMD=-(COS(D)*D-SIN(D))/D**3*.5 ELSE POM=1. POMD=0. ENDIF IF(M.EQ.0) THEN F=D1*(POM-1.)+1. DF1=D1*POMD DF2=POM-1. ELSE ARG=PI2*FLOAT(M)*.5 ALFA1=ARG*D1 CS1=COS(ALFA1) SN1=SIN(ALFA1) FCIT=SN1/ARG F=(POM-1.)*FCIT DF1=POMD*FCIT DF2=(POM-1.)*CS1 DF3=0. ENDIF RETURN 3000 IF(M.EQ.0) THEN F=D1 DF1=0. DF2=1. DF3=0. ELSE ARG=PI2*FLOAT(M)*.5 ALFA1=ARG*D1 CS1=COS(ALFA1) SN1=SIN(ALFA1) F=SN1/ARG DF1=0. DF2=CS1 DF3=0. ENDIF RETURN END SUBROUTINE XMOD(F,DF1,DF2,DF3,U,U1,U2,M,ITYP) COMMON/STALY/ PI2 POM=PI2*FLOAT(M)*.5*U1 ARG=POM+U IF(ARG.NE.0.) THEN ARG2=ARG**2 SN=SIN(ARG) CS=COS(ARG) F=U1*SN/ARG DF1=U1*(CS*ARG-SN)/ARG2 DF2=((SN+POM*CS)*ARG-POM*SN)/ARG2 ELSE F=U1 DF1=0. DF2=1. ENDIF DF3=0. RETURN END /EXEC $TSOSLNK PROG PXXX,FILENAM=P.JANA4S.POM,IDA=N INCLUDE * RESOLVE ,$X.USER.FOR1 END /COPY P.JANA4S.POM,P.JANA4S /ER P.JANA4S.POM /SYSFILE SYSLST=(PRIMARY) /EXEC $X.P.SRCXREF C.JANA4S.POM C.JANA4S END /ER C.JANA4S /CAT C.JANA4S.POM,C.JANA4S,STATE=U /LOGOFF NOSPOOL