C-----PROFILE REFINEMENT PROGRAM C-----MODIFIED DECEMBER 1972 FOR ANISOTROPIC TEMPERATURE FACTORS BY C-----A.W.HEWAT, UKAERE HARWELL, DIDCOT, BERKSHIRE, ENGLAND C-----THE SUBPROGRAMS CALCUL,SUMMAT AND FLEES ARE CALLED MOST FREQUENTLY C-----AND SHOULD THEREFORE BE KEPT IN CORE PERMANENTLY. C-----H.M.RIETVELD-DECEMBER 1970. C ******************************************************************* C-----VAX-8650, 21 Jan 1987. ISS,IALFA,IOMEG,ISTAP --> YSS,YALFA,YOMEGA,YSTAP C C*****INCLUDES VAX FLOATING OVERFLOW TRAP (EXTERNAL HAND_DIV0) C***************A.W.HEWAT 24-JUN-86************************************* C C * VARIABLE DIMENSIONS * C C * R=MAXIMUM NUMBER OF EQUIVALENT POSITIONS * C * K=MAXIMUM NUMBER OF SCATTERING LENGTHS * C * F=MAXIMUM NUMBER OF FORMFACTORS * C * N=MAXIMUM NUMBER OF ATOMS * C * M=MAXIMUM NUMBER OF MAGNETIC ROTATION MATRICES * C * P=MAXIMUM NUMBER OF LEAST SQUARES PARAMETERS * C * C=MAXIMUM NUMBER OF CONSTRAINT FUNCTIONS * C * T=MAXIMUM NUMBER OF TERMS IN CONSTRAINT FUNCTION * C * L=MAXIMUM NUMBER OF OVERLAPPING REFLECTIONS * C * I=MAXIMUM NUMBER OF REFLECTIONS * C * R1=R+1 * C * N2=N+2 * C * PC=P+C * C * ***** * C * ARRAYS WITH VARIABLE DIMENSIONS * C * ICODE(I) MEQ(N) MTYP(N) * C * NTYP(N) ICON(C,T) JCON(C,T) * C * LCONI(C,T) LCONJ(C,T) LN(C) * C * LP(N2,14) A(N2,14) ATEXT(N) * C * B(K) C(C,L) CF(C) * C * CONA(C,T) COSAR(N,R1) COSARI(R1) * C * DERIV(PC) DERSTO(L,PC) DF(C) * C * DUMP(L) F(F,21) FAC(I) * C * FESD(I) FF(I,F) FLAG(N2,14) * C * FMAG(I) FNUC(I) FOBS(I) * C * H(R1,3) HH(I,3) HW(L) * C * OCC(N) PAK(I) RM(R1,M,9) * C * RSA(N,3,3) RSB(N,3,3) S(F,21) * C * SA(N) SB(N) SINAR(N,R1) * C * SINARI(R1) SJJ(I) SM(R1,3,3) * C * SMM(PC,PC) SNEX(N) SOME(L) * C * SOMEGA(I) T(R1,3) TANN(L) * C * TEMP(N) THETA(I) TR(R1) * C * V(PC) XL(N2,14) SZ(N2,14) * C ******************************************************************* C * C ****************** C * PRESENT VALUES * C * R =33 * C * K = 6 * C * F = 2 * C * N= 69 * C * N2=71 * C * M = 9 * C * P =145 * C * C = 4 * C * PC=149 * C * L =98 * C * I =4889 * C * T = 9 * C * R1=34 * C ****************** C **************************** C SUBROUTINE STRUCTURE C MAIN 930 - INPUT 480 - FLEES -I C 1620 - CELL I C 2570 - PROFIL -I I C MAIN 990 - ITER 440 - CONPRE I I C 560 - FLEES I-I C 660 - PROFIL -I C 930 - CALCUL C 950 - SUMMAT C MAIN1000 - MATRIX370 - DTSLIN C 1010 - DIRECT840 - ESD 260 - ERROR C 1020 - OUTPUT C 1040 - EXPUT C ****************************************************************************** C CHANGES FOR RHEL SYSTEM C ****************************************************************************** C * UNIT 6 = LINEPRINTER * C * UNIT 5 = CARDREADER * C * UNIT 7 = CARDPUNCH * C * UNIT 2 = TEMPORARY DATA SET ON DISC * C **************************** LOGICAL HALPER,PAC,UNIAX,CUBIC,PUNCH,TYPE,STAR,SKIP, 1YES,VERT,FORCED COMMON CUBIC,FORCED,HALPER,PAC,PUNCH,SKIP,STAR,TYPE,UNIAX, 1VERT,YES COMMON AH,B1,B2,BH,CH,DELTA,EPS,FACCOS,FACCUB,FACSIN,FLOG,FLOG2, 1FMAGDI,FMAGOB,FNUCDI,FNUCOB,OMEGA,PI,RELAXB,RELAXC,RELAXS,SCALE, 2SKEW,SOM,SUMCAL,SUMDFF,SUMDIF,SUMFOB,SUMOBS,SUMW,SUMWBS,SW,W,XN, 3YCALC,YOBS,ZERO,SLABDA,SLOG,ARG COMMON YALFA,ICENT,ICYCLE,YOMEG,IORD1,IORD2,IPONS,IPREV,IRL,YSTAP, 1KL,KM,KSTART,LCOORD,LCORD,LIM,LIMO,LLL,LMAT,LOUT,MAXCON,MAXS,MN,N, 2NCYCLE,NKEY,NLC,NN,NQC,NTEXT,NUM,YSS,NO COMMON ICODE(4889),ICON(4,9),JCON(4,9),LCONI(4,9),LCONJ(4,9),LN(4) 1,LP(71,14),MEQ(69),MTYP(69),NTYP(69),IOBS(10),ICALC(10) COMMON A(71,14),AL(3,3),ATEXT(69),B(6),C(4,98),CF(4),CONA(4,9), 1COSA(3),COSAF(3),COSAR(69,48),COSARI(48),DERIV(149),DERSTO(98,149) 2,DF(4),DUMP(98),DUMPA(500),F(2,21),FAC(4889),FESD(4889),FF(4889,2) 3,FLAG(71,14),FMAG(4889),FNUC(4889),FOBS(4889),H(48,3),HH(4889,3), 4HL(3),HN(3),HW(98),OCC(69),PAK(4889),PH(3),RM(48,9,9),RSA(69,3,3), 5RSB(69,3,3),S(2,21),SA(69),SB(69),SIDE(3),SINAR(69,48),SINARI(48), 6SJJ(4889),SM(48,3,3),SMM(149,149),SMOM(3),SNEX(69),SOME(98), 7SOMEGA(4889),T(48,3),TANN(98),TEMP(69),TEXT(20),THETA(4889),TR(48) 8,V(149),XL(71,14),SZ(71,14) CHARACTER*20 FILEN INTEGER*4 CLOCK_TICKS,CLOCK_TICKS2,SECONDS COMMON NRD,NWD,NSD,NID,NFD,NDD,NCD,PHASEN(4889),SP(6,6),RFAC(6) 1,URABS,ABSOR COMMON/CPUTIM/START_TIME,CLOCK_TICKS C **** Assume all variables initialised to zero DO 1001I=1,PC 1001 V(I)=0.0 ! A.Hewat 1/2/89 otherwise problems with OUTPUT DO65J=1,4 etc NRD=30 NWD=31 NSD=2 NID=33 NFD=34 NDD=38 NCD=39 C DATA NRD,NWD,NSD,NID,NFD,NDD,NCD/30,31,2,33,34,38,39/ WRITE(*,1009) 1009 FORMAT(' PROF2 Rietveld refinement. **PROFF contains F not F**2', A //' File name for profile refinement ? eg PROFD.DAT ') READ(*,1019) FILEN 1019 FORMAT(A) IF(FILEN.EQ.' ') FILEN='PROFD.DAT' C.... Open statements modified. J.K.C. OPEN(UNIT=NSD,FILE='PROF1.BIN',STATUS='OLD',FORM='UNFORMATTED') OPEN(UNIT=NRD,FILE=FILEN,STATUS='OLD',SHARED,READONLY) OPEN(UNIT=NWD,FILE='PROF2.LIS',STATUS='NEW') OPEN(UNIT=NID,FILE='PROFI.DAT',STATUS='NEW', + CARRIAGE CONTROL='LIST') OPEN(UNIT=NFD,FILE='PROFF.DAT',STATUS='NEW', + CARRIAGE CONTROL='LIST') OPEN(UNIT=NDD,FILE='PROFD.DAT',STATUS='NEW', + CARRIAGE CONTROL='LIST') OPEN(UNIT=NCD,FILE='PROFC.DAT',STATUS='NEW', + CARRIAGE CONTROL='LIST') WRITE(*,1029) 1029 FORMAT(' Please wait...Profile refinement in progress') c C cpu timer *** A.Hewat 23 Feb 89 *** SECONDS=0 CALL TIMER(SECONDS) IF(SECONDS.NE.0) START_TIME=SECONDS C cpu timer *** A.Hewat 23 Feb 89 *** C C-----WHEN ANGLES IN 1/10000 CYCLES, FOLLOWING STATEMENT SHOULD READ: C-----PI=10000./3.14159265359 PI=36000./3.14159265359 1000 CALL INPUT IF(NKEY.EQ.12345)GOTO 3 ICYCLE=0 FORCED=.FALSE. 1 ICYCLE=ICYCLE+1 IF(ICYCLE.GT.NCYCLE)GOTO 2 CALL ITER IF(MAXS.NE.0)CALL MATRIX IF(TYPE)CALL DIRECT CALL OUTPUT GOTO 1 2 CALL EXPUT 3 CONTINUE CLOSE(UNIT=NSD) CLOSE(UNIT=NRD) CLOSE(UNIT=NWD) CLOSE(UNIT=NID) CLOSE(UNIT=NFD) CLOSE(UNIT=NDD) CLOSE(UNIT=NCD) C cpu timer *** A.Hewat 23 Feb 89 *** SECONDS=1 CALL TIMER(SECONDS) IF(SECONDS.EQ.1) GO TO 1234 ETIME=FLOAT(SECONDS)-START_TIME CTIME=ETIME Write(*,9999) CTIME,ETIME 9999 FORMAT(//' CPU time :',F9.2,5X,'Elapsed time :',F9.2/) 1234 PAUSE ' Press RETURN to finish.' C cpu timer *** A.Hewat 23 Feb 89 *** STOP END SUBROUTINE INPUT C-----THIS SUBROUTINE READS AND PRINTS THE INPUT DATA,AND INITIALIZES C-----OTHER DATA COMMON CUBIC,FORCED,HALPER,PAC,PUNCH,SKIP,STAR,TYPE,UNIAX, 1VERT,YES COMMON AH,B1,B2,BH,CH,DELTA,EPS,FACCOS,FACCUB,FACSIN,FLOG,FLOG2, 1FMAGDI,FMAGOB,FNUCDI,FNUCOB,OMEGA,PI,RELAXB,RELAXC,RELAXS,SCALE, 2SKEW,SOM,SUMCAL,SUMDFF,SUMDIF,SUMFOB,SUMOBS,SUMW,SUMWBS,SW,W,XN, 3YCALC,YOBS,ZERO,SLABDA,SLOG,ARG COMMON YALFA,ICENT,ICYCLE,YOMEG,IORD1,IORD2,IPONS,IPREV,IRL,YSTAP, 1KL,KM,KSTART,LCOORD,LCORD,LIM,LIMO,LLL,LMAT,LOUT,MAXCON,MAXS,MN,N, 2NCYCLE,NKEY,NLC,NN,NQC,NTEXT,NUM,YSS,NO COMMON ICODE(4889),ICON(4,9),JCON(4,9),LCONI(4,9),LCONJ(4,9),LN(4) 1,LP(71,14),MEQ(69),MTYP(69),NTYP(69),IOBS(10),ICALC(10) COMMON A(71,14),AL(3,3),ATEXT(69),B(6),C(4,98),CF(4),CONA(4,9), 1COSA(3),COSAF(3),COSAR(69,48),COSARI(48),DERIV(149),DERSTO(98,149) 2,DF(4),DUMP(98),DUMPA(500),F(2,21),FAC(4889),FESD(4889),FF(4889,2) 3,FLAG(71,14),FMAG(4889),FNUC(4889),FOBS(4889),H(48,3),HH(4889,3), 4HL(3),HN(3),HW(98),OCC(69),PAK(4889),PH(3),RM(48,9,9),RSA(69,3,3), 5RSB(69,3,3),S(2,21),SA(69),SB(69),SIDE(3),SINAR(69,48),SINARI(48), 6SJJ(4889),SM(48,3,3),SMM(149,149),SMOM(3),SNEX(69),SOME(98), 7SOMEGA(4889),T(48,3),TANN(98),TEMP(69),TEXT(20),THETA(4889),TR(48) 8,V(149),XL(71,14),SZ(71,14) COMMON NRD,NWD,NSD,NID,NFD,NDD,NCD,PHASEN(4889),SP(6,6),RFAC(6) 1,URABS,ABSOR LOGICAL HALPER,PAC,UNIAX,CUBIC,PUNCH,TYPE,STAR,SKIP, 1YES,VERT,FORCED LLL=0 C-----READ PROBLEM IDENTIFICATION READ(NRD,1)TEXT C-----READ PROBLEM PARAMETERS READ(NRD,2)FILEN,EPS,K,(PH(I),I=1,3),LIM IF((ABS(PH(1))+ABS(PH(2))+ABS(PH(3))).LE.0.001) PH(3)=1. !A.Hewat 31/1/89 MAGTYP=K READ(NRD,3)ICENT,IRL,KL,KM,N,MN,SLABDA IRL=IRL+1 C-----PRINT PROGRAM HEADING WRITE(NWD,5) IF(LIM.LT.0) WRITE(NWD,401) CJH WRITE(NWD,160)TEXT C-----SEARCH FOR FILE NAME C ****************************************************************************** C CHANGES FOR RHEL SYSTEM, ONLY ONE DATA SET ON DISC 213 READ(NSD,END=212)DUMPA IF(DUMPA(2).EQ.FILEN)GOTO 211 C END OF THIS SET OF CHANGES C ****************************************************************************** GOTO 213 211 WRITE(NWD,9)FILEN LLL=2 NKEY=1 GOTO 214 212 WRITE(NWD,215)FILEN X=FILEN REWIND 2 NKEY=12345 GOTO 221 214 LIMO=FLEES(0.)+0.5 WRITE(NWD,10)EPS HALPER=.FALSE. UNIAX=.FALSE. CUBIC=.FALSE. IF(K-1)11,12,11 12 WRITE(NWD,13) HALPER=.TRUE. 11 IF(K-2)14,15,14 15 WRITE(NWD,16) UNIAX=.TRUE. HL(1)=0. HL(2)=0. HL(3)=1. 14 IF(K-3)17,18,17 18 WRITE(NWD,19) CUBIC=.TRUE. FACCUB=2./3. 17 WRITE(NWD,20)(PH(I),I=1,3) WRITE(NWD,21)SLABDA SLAMDA=SLABDA SLABDA=0.25*SLABDA*SLABDA PI=SLABDA*PI C-----READ EQUIVALENT POSITION VECTOR AND MATRICES DO 22 I=1,3 T(1,I)=0. DO 22 J=1,3 IJ=(I-1)*3+J IF(I-J)23,24,23 24 SM(1,I,J)=1. GOTO 25 23 SM(1,I,J)=0. IF(MN)22,22,25 25 DO 180 K=1,MN IF(I-J)26,27,26 27 RM(1,K,IJ)=1. GOTO 180 26 RM(1,K,IJ)=0. 180 CONTINUE 22 CONTINUE IF(IRL-1)229,229,28 28 DO 29 IR=2,IRL READ(NRD,30)((SM(IR,I,J),J=1,3),T(IR,I),I=1,3) IF(MN)29,29,31 31 DO 181 K=1,MN 181 READ(NRD,30)(RM(IR,K,IJ),IJ=1,9) 29 CONTINUE C-----READ SCATTERING LENGTHS 229 READ(NRD,33)(B(I),I=1,KL) IF(KM)34,34,35 C-----READ FORMFACTORS 35 DO 36 I=1,KM TYPE=.TRUE. DO 36 J=1,21 IF(.NOT.TYPE)GOTO 38 NFF=J-1 READ(NRD,37)S(I,J),F(I,J) IF(S(I,J)+100.)36,39,36 39 TYPE=.FALSE. S(I,J)=S(I,J-1) X=S(I,J) F(I,J)=F(I,J-1) Y=F(I,J) GOTO 36 38 S(I,J)=X F(I,J)=Y 36 CONTINUE C-----READ NUMBER OF CYCLES,RELAXATION FACTORS AND OUTPUT PARAMETERS 34 READ(NRD,40) A NCYCLE,RELAXC,RELAXB,RELAXS,LOUT,IPONS,LMAT,LCORD,LCOORD C IF REQUIRED, PREPARE NEW COORDINATES FILE IF(LCOORD) 348,348,340 340 WRITE(NDD,1) TEXT WRITE(NDD,2) FILEN,EPS,MAGTYP,(PH(I),I=1,3),LIM I=IRL-1 WRITE(NDD,3) ICENT,I,KL,KM,N,MN,SLAMDA IF(IRL-1) 345,345,341 341 FLAGD=' ' DO344IR=2,IRL IF(MN.EQ.0.AND.IR.EQ.IRL) FLAGD='/' WRITE(NDD,30)((SM(IR,I,J),J=1,3),T(IR,I),I=1,3),FLAGD IF(MN) 344,344,342 342 DO343K=1,MN IF(K.EQ.MN.AND.IR.EQ.IRL) FLAGD='/' 343 WRITE(NDD,319) (RM(IR,K,IJ),IJ=1,9),FLAGD 344 CONTINUE 345 WRITE(NDD,33) (B(I),I=1,KL) IF(KM) 348,348,346 346 DO3462 I=1,KM DO3461 J=1,20 WRITE (NDD,37)S(I,J),F(I,J) IF((S(I,J).EQ.S(I,J+1)).AND.(F(I,J).EQ.F(I,J+1)))GOTO3462 3461 CONTINUE 3462 WRITE(NDD,349) 349 FORMAT(8H-100. ,71X,1H/) C-----PRINT ALL INPUT DATA 348 WRITE(NWD,41)NCYCLE IF(LOUT.EQ.1)WRITE(NWD,42) IF(LMAT.EQ.1)WRITE(NWD,43) IF((IPONS.EQ.1).OR.(IPONS.EQ.3))WRITE(NWD,44) IF((IPONS.EQ.2).OR.(IPONS.EQ.3))WRITE(NWD,45) IF(LCOORD.EQ.1)WRITE(NWD,46) IF(LCORD.EQ.1)WRITE(NWD,47) WRITE(NWD,48)RELAXC,RELAXB,RELAXS IF(ICENT.EQ.2)GOTO 300 WRITE(NWD,301) GOTO 302 300 WRITE(NWD,303) 302 WRITE(NWD,49) IF(IRL-2)53,51,51 51 DO 52 IR=2,IRL WRITE(NWD,152)((SM(IR,I,J),J=1,3),T(IR,I),I=1,3) IF(MN)52,52,54 54 DO182 K=1,MN 182 WRITE(NWD,55)(RM(IR,K,IJ),IJ=1,9) 52 CONTINUE 53 WRITE(NWD,50) WRITE(NWD,56)(I,B(I),I=1,KL) IF(KM)260,260,58 58 WRITE(NWD,59) DO 60 I=1,KM WRITE(NWD,61) DO 63 J=1,20 WRITE(NWD,62)S(I,J),F(I,J) IF((S(I,J).EQ.S(I,J+1)).AND.(F(I,J).EQ.F(I,J+1)))GOTO 60 63 CONTINUE 60 CONTINUE 260 WRITE(NWD,57) C-----READ SCALE AND TEMP FACTOR and absorbtion uR !A.Hewat 4 May'90 READ(NRD,64)SCALE,XL(N+1,4),URABS !A.Hewat 4 May'90 XL(N+1,5)=1./SCALE ABSOR=-0.5*SLAMDA*SLAMDA*(0.0368+0.3750*URABS)*URABS C-----READ AND PRINT FOR EACH ATOM READ(NRD,65)(ATEXT(I),NTYP(I),MTYP(I),MEQ(I), A (XL(I,J),J=1,14),I=1,N) WRITE(NWD,170)(ATEXT(I),NTYP(I),MTYP(I),MEQ(I),(XL(I,J),J=1,14), 1I=1,N) WRITE(NWD,66)SCALE,XL(N+1,4),URABS,ABSOR C-----READ AND PRINT PROFILE PARAMETERS READ(NRD,67)(XL(N+2,I),I=1,4) WRITE(NWD,68)(XL(N+2,I),I=1,3) WRITE(NWD,69)XL(N+2,4) READ(NRD,70)I IF(I)71,72,72 71 CALL CELL DO 73 I=1,3 73 XL(N+2,I+4)=AL(I,I) XL(N+2,8)=AL(2,3) XL(N+2,9)=AL(1,3) XL(N+2,10)=AL(1,2) GOTO 74 72 READ(NRD,75)(XL(N+2,I),I=5,10) 74 WRITE(NWD,76)(XL(N+2,I),I=5,10) READ(NRD,77)XL(N+1,1),XL(N+1,2) WRITE(NWD,78)(XL(N+1,I),I=1,2) N2=N+2 DO 79 I=1,N2 DO 79 J=1,14 A(I,J)=0. FLAG(I,J)=0. 79 LP(I,J)=0 MAXCON=0 READ(NRD,70)MAXS IF(MAXS)80,81,80 80 WRITE(NWD,82) C-----READ AND PRINT CODEWORDS FOR ATOMIC PARAMETERS DO 83 I=1,N READ(NRD,86)(FLAG(I,J),J=1,14) DO 84 J=1,14 X=FLAG(I,J) Y=INT(ABS(X))/10 LP(I,J)=Y 84 A(I,J)=(ABS(X)-10.*Y)*SIGN(1.,X) WRITE(NWD,85)ATEXT(I),(FLAG(I,J),J=1,14) 83 CONTINUE C-----READ AND PRINT CODEWORDS FOR PROFILE PARAMETERS READ(NRD,86)FLAG(N+1,5),FLAG(N+1,4) READ(NRD,86)(FLAG(N+2,J),J=1,4) READ(NRD,86)(FLAG(N+2,J),J=5,10) READ(NRD,86)(FLAG(N+1,J),J=1,2) DO 87 I=1,5 X=FLAG(N+1,I) Y=INT(ABS(X))/10 LP(N+1,I)=Y 87 A(N+1,I)=(ABS(X)-10.*Y)*SIGN(1.,X) DO 88 I=1,10 X=FLAG(N+2,I) Y=INT(ABS(X))/10 LP(N+2,I)=Y 88 A(N+2,I)=(ABS(X)-10.*Y)*SIGN(1.,X) WRITE(NWD,89)FLAG(N+1,5) WRITE(NWD,90)FLAG(N+1,4) WRITE(NWD,91)(FLAG(N+2,J),J=1,3) WRITE(NWD,92)FLAG(N+2,4) WRITE(NWD,93)(FLAG(N+2,J),J=5,10) WRITE(NWD,95)FLAG(N+1,1) WRITE(NWD,94)FLAG(N+1,2) C-----READ CONSTRAINTS READ(NRD,96)NLC,NQC MAXCON=NLC+NQC IF(MAXCON.EQ.0)GOTO 81 IL=0 DO 97 I=1,MAXCON 103 IL=IL+1 READ(NRD,98) A CONA(I,IL),LCONI(I,IL),LCONJ(I,IL),ICON(I,IL),JCON(I,IL) IF(CONA(I,IL)+100.)99,100,99 99 IF(I-NLC)101,101,102 101 ICON(I,IL)=LCONI(I,IL) JCON(I,IL)=LCONJ(I,IL) 102 GOTO 103 100 LN(I)=IL-1 READ(NRD,98)DF(I) IL=0 97 CONTINUE 81 IF(MAXCON.NE.0)GOTO 201 WRITE(NWD,202) GOTO 200 201 WRITE(NWD,203)MAXCON IF(NLC.EQ.0)GOTO 204 DO 205 I=1,NLC LNI=LN(I) WRITE(NWD,207)(CONA(I,J),ICON(I,J),JCON(I,J),J=1,LNI) 205 WRITE(NWD,208) DF(I) IF(NQC.EQ.0)GOTO 200 204 DO 209 IL=1,NQC I=IL+NLC WRITE(NWD,206) LNI=LN(I) WRITE(NWD,210) A (CONA(I,J),LCONI(I,J),LCONJ(I,J),ICON(I,J),JCON(I,J), 1J=1,LNI) 209 WRITE(NWD,208)DF(I) 200 PAC=.NOT.((XL(N+1,1).EQ.0.).AND.(FLAG(N+1,1).EQ.0.)) C-----READ REFLECTIONS FROM MAGNETIC TAPE NN=FLEES(0.)+0.5 DO 104 I=1,3 104 HN(I)=FLEES(0.) FLOG=4.*ALOG(2.) FLOG2=2.*FLOG SLOG=1./SQRT(FLOG) CALL PROFIL IF(.NOT.PAC)GOTO 105 C-----CALCULATE FOR EACH REFLECTION (1+COS(ETA)**2)/2 DO 106 I=1,3 106 PH(I)=PH(I)/HN(I) TT=0. DO 107 I=1,3 DO 107 J=I,3 107 TT=AL(I,J)*PH(I)*PH(J)+TT 105 IF(.NOT.UNIAX)GOTO 108 C-----CALCULATE FOR EACH REFLECTION ALPHA**2 DO 109 I=1,3 109 HL(I)=HL(I)/HN(I) Z=0. DO 110 I=1,3 DO 110 J=I,3 110 Z=AL(I,J)*HL(I)*HL(J)+Z 108 DO 121 I=1,NN ICODE(I)=FLEES(0.)+0.5 DO 112 J=1,3 112 HH(I,J)=FLEES(0.)/HN(J) SJJ(I)=FLEES(0.) X=0. DO 113 J=1,3 DO 113 K=J,3 113 X=AL(J,K)*HH(I,J)*HH(I,K)+X IF(.NOT.(UNIAX.AND.(ICODE(I).NE.1)))GOTO 114 Y=0. DO 115 J=1,3 DO 115 K=J,3 115 Y=AL(J,K)*(HL(J)*HH(I,K)+HL(K)*HH(I,J))*0.5+Y FAC(I)=Y*Y/(Z*X) 114 IF(.NOT.PAC)GOTO 116 BB=0. DO 117 J=1,3 DO 117 K=J,3 117 BB=AL(J,K)*(PH(J)*HH(I,K)+PH(K)*HH(I,J))*0.5+BB BB=BB*BB/(TT*X) IF(BB)118,119,118 119 BB=1.5707963268 GOTO 120 118 BB=ATAN(SQRT(ABS((1.-BB)/BB))) 120 PAK(I)=BB*BB 116 X=SQRT(X)*0.5 IF(KM)121,121,122 122 DO 183 J=1,KM DO 123 K=1,20 K1=K-1 IF(K1.EQ.0)K1=1 IF(X.GT.S(J,K))GOTO 123 FF(I,J)=(S(J,K1)-X)/(S(J,K1)-S(J,K))*(F(J,K)-F(J,K1))+F(J,K1) GOTO 183 123 CONTINUE FF(I,J)=F(J,20) 183 CONTINUE 121 CONTINUE 221 CONTINUE 1 FORMAT(20A4) 2 FORMAT(A4,4X,F8.1,I8,3F8.0,I8,23X,1H/) 3 FORMAT(6I8,F8.4,23X,1H/) 5 FORMAT(' NEUTRON POWDER PROFILE REFINEMENT FOR NUCLEAR & MAGNETIC 1STRUCTURES -ILL VERSION 149 PARAMETERS & 4889 REFLEXIONS FEB 78'/ 2 ' ***** ALL PUBLISHED USE SHOULD REFER TO...'/ 3 ' 1. RIETVELD H.M (1969) J APPL.CRYST 2,65'/ 4 ' 2. HEWAT A.W.(1973) HARWELL REPORT AERE-R7350'/) 401 FORMAT('+ 3. HOWARD C.J. (1981) J APPL.CRYST. 15,615'/) CJH 9 FORMAT(12H FILE NAMED ,A4,1X,8HLOCATED.) 10 FORMAT(11H EPS-VALUE=,F6.1) 13 FORMAT(67H MAGNETIC INTENSITIES CALCULATION ACCORDING TO HALPERN AHHHHHHHH 1ND JOHNSON.) 16 FORMAT(41H UNIAXIAL CONFIGURATIONAL SPIN SYMMETRY.) 19 FORMAT(37H CUBIC CONFIGURATIONAL SPIN SYMMETRY.) 20 FORMAT(33H PREFERRED ORIENTATION DIRECTION=,3F3.0) 21 FORMAT(12H WAVELENGTH=,F8.4) 30 FORMAT(3(4F6.2),7X,A1) 319 FORMAT(3(3F6.2),25X,A1) 32 FORMAT(9F8.2) 33 FORMAT(6F8.3,31X,1H/) 37 FORMAT(2F8.4) 40 FORMAT(I8,3F8.2,5I8) 41 FORMAT(18H NUMBER OF CYCLES=,I6) 42 FORMAT(55H PRINTED OUTPUT OF OBS + CALC INTENSITIES ON LAST CYCLE) 43 FORMAT(55H PRINTED OUTPUT OF CORRELATION - MATRIX ON LAST CYCLE) 44 FORMAT(55H PUNCHED OUTPUT OF STRUCTURE FACTOR TAPE ON LAST CYCLE) 45 FORMAT(55H PUNCHED OUTPUT OF OBS + CALC INTENSITIES ON LAST CYCLE) 46 FORMAT(55H PUNCHED OUTPUT OF NEW COORDINATES TAPE ON LAST CYCLE) 47 FORMAT(55H PUNCHED OUTPUT COORDINATES COVARIANCES ON LAST CYCLE) 48 FORMAT(19H0RELAXATION FACTORS,/21H FOR COORDINATES =,F5.2,/21H HHHHHHHH 1FOR TEMP FACTORS =,F5.2,/21H FOR SCALING FACTORS=,F5.2) 49 FORMAT(22H0EQUIVALENT POSITIONS=) 50 FORMAT(20H SCATTERING LENGTHS=) 55 FORMAT(5X,3F5.1,/5H M= ,3F5.1,/,5X,3F5.1/) 56 FORMAT(3H B(,I1,2H)=,F8.3) 57 FORMAT(/25H ***INITIAL PARAMETERS***/18H ATOM NTP MTP MRT,3X,1HX,INPU3450 17X,1HY,7X,1HZ,7X,1HB,7X,1HN,6X,2HKX,6X,2HKY,6X,2HKZ,5X,3HB11,5X, INPU3460 23HB22,5X,3HB33,5X,3HB12,5X,3HB13,5X,3HB23) INPU3461 59 FORMAT(13H0FORMFACTORS=) INPU3470 61 FORMAT(20H SINTH/LABDA F) INPU3480 62 FORMAT(1X,F8.4,F13.4) INPU3490 64 FORMAT(3F8.6) INPU3500 65 FORMAT(A4,3I4,8F8.5/6F8.5) INPU3510 66 FORMAT(23H OVERALL SCALE FACTOR =,F12.5/23H OVERALL TEMP. FACTOR =INPU3520 1,F12.5/23H OVERALL ABSORBTION UR=,F12.5/23H DELTA-B ABSORBTION = 2,F12.5) INPU3530 67 FORMAT(4F8.2) INPU3540 68 FORMAT(22H HALFWIDTH PARAMETERS=,3F12.4) INPU3550 69 FORMAT(11H ZEROPOINT=,F7.3) INPU3560 70 FORMAT(I8) INPU3570 75 FORMAT(6F8.4) INPU3580 76 FORMAT(16H CELL CONSTANTS=,6F12.7) INPU3590 77 FORMAT(2F8.4) INPU3600 78 FORMAT(33H PREFERRED ORIENTATION PARAMETER=,F7.3,/,21H ASYMMETRY PINPU3610 1ARAMETER=,F10.6) INPU3620 82 FORMAT(26H0***CODING OF VARIABLES***/5H ATOM,16X,1HX,7X,1HY,7X,1HZINPU3630 1,7X,1HB,7X,1HN,6X,2HKX,6X,2HKY,6X,2HKZ,5X,3HB11,5X,3HB22,5X,3HB33,INPU3640 25X,3HB12,5X,3HB13,5X,3HB23) 1 85 FORMAT(1X,A4,10X,14F8.3) INPU3650 86 FORMAT(8F8.2) INPU3660 89 FORMAT(22H OVERALL SCALE FACTOR=,F8.2) INPU3670 90 FORMAT(22H OVERALL TEMP. FACTOR=,F8.2) INPU3680 91 FORMAT(22H HALFWIDTH PARAMETERS=,3F8.2) INPU3690 92 FORMAT(11H ZEROPOINT=,F8.2) INPU3700 93 FORMAT(16H CELL CONSTANTS=,6F8.2) INPU3710 94 FORMAT(21H ASYMMETRY PARAMETER=,F8.2) INPU3720 95 FORMAT(33H PREFERRED ORIENTATION PARAMETER=,F8.2) INPU3730 96 FORMAT(2I8) INPU3740 98 FORMAT(F8.4,4I8) INPU3750 152 FORMAT(5X,3F5.1,8X,F5.2,/,5H R= ,3F5.1,8H T= ,F5.2,/,5X,3F5.1INPU3760 1,8X,F5.2/) INPU3770 160 FORMAT(/4H ***,20A4,3H***/) INPU3780 170 FORMAT(1X,A4,3I4,14F8.5) INPU3790 202 FORMAT(25H NO CONSTRAINT FUNCTIONS.) INPU3800 203 FORMAT(1X,I2,1X,21HCONSTRAINT FUNCTIONS=) INPU3810 206 FORMAT(/) INPU3820 207 FORMAT(8(F6.2,3H*X(,I2,1H,,I2,1H))) INPU3830 208 FORMAT(3H D=,F6.2) INPU3840 210 FORMAT(5(F6.2,3H*X(,I2,1H,,I2,5H)*DX(,I2,1H,,I2,1H))) INPU3850 215 FORMAT(12H FILE NAMED ,A4,1X,12HNOT LOCATED.) INPU3860 301 FORMAT(38H0THE STRUCTURE IS NON CENTROSYMMETRIC.) INPU3870 303 FORMAT(34H0THE STRUCTURE IS CENTROSYMMETRIC.) INPU3880 RETURN INPU3890 END INPU3900 FUNCTION FLEES(AA) C-----THIS FUNCTION CAUSES A TRANSFER OF 500 WORDS FROM TAPE INTO THE C-----ASSIGNING CONSECUTIVE ELEMENT VALUES TO THE FUNCTION. COMMON CUBIC,FORCED,HALPER,PAC,PUNCH,SKIP,STAR,TYPE,UNIAX, 1VERT,YES COMMON AH,B1,B2,BH,CH,DELTA,EPS,FACCOS,FACCUB,FACSIN,FLOG,FLOG2, 1FMAGDI,FMAGOB,FNUCDI,FNUCOB,OMEGA,PI,RELAXB,RELAXC,RELAXS,SCALE, 2SKEW,SOM,SUMCAL,SUMDFF,SUMDIF,SUMFOB,SUMOBS,SUMW,SUMWBS,SW,W,XN, 3YCALC,YOBS,ZERO,SLABDA,SLOG,ARG COMMON YALFA,ICENT,ICYCLE,YOMEG,IORD1,IORD2,IPONS,IPREV,IRL,YSTAP, 1KL,KM,KSTART,LCOORD,LCORD,LIM,LIMO,LLL,LMAT,LOUT,MAXCON,MAXS,MN,N, 2NCYCLE,NKEY,NLC,NN,NQC,NTEXT,NUM,YSS,NO COMMON ICODE(4889),ICON(4,9),JCON(4,9),LCONI(4,9),LCONJ(4,9),LN(4) 1,LP(71,14),MEQ(69),MTYP(69),NTYP(69),IOBS(10),ICALC(10) COMMON A(71,14),AL(3,3),ATEXT(69),B(6),C(4,98),CF(4),CONA(4,9), 1COSA(3),COSAF(3),COSAR(69,48),COSARI(48),DERIV(149),DERSTO(98,149) 2,DF(4),DUMP(98),DUMPA(500),F(2,21),FAC(4889),FESD(4889),FF(4889,2) 3,FLAG(71,14),FMAG(4889),FNUC(4889),FOBS(4889),H(48,3),HH(4889,3), 4HL(3),HN(3),HW(98),OCC(69),PAK(4889),PH(3),RM(48,9,9),RSA(69,3,3), 5RSB(69,3,3),S(2,21),SA(69),SB(69),SIDE(3),SINAR(69,48),SINARI(48), 6SJJ(4889),SM(48,3,3),SMM(149,149),SMOM(3),SNEX(69),SOME(98), 7SOMEGA(4889),T(48,3),TANN(98),TEMP(69),TEXT(20),THETA(4889),TR(48) 8,V(149),XL(71,14),SZ(71,14) COMMON NRD,NWD,NSD,NID,NFD,NDD,NCD,PHASEN(4889),SP(6,6),RFAC(6) 1,URABS,ABSOR LOGICAL HALPER,PAC,UNIAX,CUBIC,PUNCH,TYPE,STAR,SKIP, 1YES,VERT,FORCED LLL=LLL+1 FLEES=DUMPA(LLL) IF(LLL.NE.500)GOTO 1 READ(NSD)DUMPA NKEY=NKEY+1 LLL=0 1 RETURN END SUBROUTINE CELL C-----THE SUBROUTINE CELL READS THE DIRECT CELL PARAMETERS AND CONVERTS C-----THEM TO THE CELL CONSTANTS,WHICH ARE PLACED IN ARRAY X COMMON CUBIC,FORCED,HALPER,PAC,PUNCH,SKIP,STAR,TYPE,UNIAX, 1VERT,YES COMMON AH,B1,B2,BH,CH,DELTA,EPS,FACCOS,FACCUB,FACSIN,FLOG,FLOG2, 1FMAGDI,FMAGOB,FNUCDI,FNUCOB,OMEGA,PI,RELAXB,RELAXC,RELAXS,SCALE, 2SKEW,SOM,SUMCAL,SUMDFF,SUMDIF,SUMFOB,SUMOBS,SUMW,SUMWBS,SW,W,XN, 3YCALC,YOBS,ZERO,SLABDA,SLOG,ARG COMMON YALFA,ICENT,ICYCLE,YOMEG,IORD1,IORD2,IPONS,IPREV,IRL,YSTAP, 1KL,KM,KSTART,LCOORD,LCORD,LIM,LIMO,LLL,LMAT,LOUT,MAXCON,MAXS,MN,N, 2NCYCLE,NKEY,NLC,NN,NQC,NTEXT,NUM,YSS,NO COMMON ICODE(4889),ICON(4,9),JCON(4,9),LCONI(4,9),LCONJ(4,9),LN(4) 1,LP(71,14),MEQ(69),MTYP(69),NTYP(69),IOBS(10),ICALC(10) COMMON Z(71,14),X (3,3),ATEXT(69),U(6),Y(4,98),CF(4),CONA(4,9), 1COSO(3),COSAF(3),COSAR(69,48),COSARI(48),DERIV(149),DERSTO(98,149) 2,DF(4),DUMP(98),DUMPA(500),G(2,21),FAC(4889),FESD(4889),FF(4889,2) 3,FLAG(71,14),FMAG(4889),FNUC(4889),FOBS(4889),H(48,3),HH(4889,3), 4HL(3),HN(3),HW(98),OCC(69),PAK(4889),PH(3),RM(48,9,9),RSA(69,3,3), 5RSB(69,3,3),S(2,21),SA(69),SB(69),SIDE(3),SINAR(69,48),SINARI(48), 6SJJ(4889),SM(48,3,3),SMM(149,149),SMOM(3),SNEX(69),SOME(98), 7SOMEGA(4889),T(48,3),TANN(98),TEMP(69),TEXT(20),THETA(4889),TR(48) 8,Q(149),XL(71,14),SZ(71,14) COMMON NRD,NWD,NSD,NID,NFD,NDD,NCD,PHASEN(4889),SP(6,6),RFAC(6) 1,URABS,ABSOR LOGICAL HALPER,PAC,UNIAX,CUBIC,PUNCH,TYPE,STAR,SKIP, 1YES,VERT,FORCED RAD=3.14159265359/180. READ(NRD,1)A,B,C,D,E,F COSA=COS(RAD*D) COSB=COS(RAD*E) COSC=COS(RAD*F) SIDE(1)=A SIDE(2)=B SIDE(3)=C COSAF(1)=COSA COSAF(2)=COSB COSAF(3)=COSC SINA=SQRT(1.-COSA*COSA) SINB=SQRT(1.-COSB*COSB) SINC=SQRT(1.-COSC*COSC) V=A*B*C*SQRT(1.-COSA*COSA-COSB*COSB-COSC*COSC+2.*COSA*COSB*COSC) AS=B*C*SINA/V BS=C*A*SINB/V CS=A*B*SINC/V COSAS=(COSB*COSC-COSA)/(SINB*SINC) COSBS=(COSC*COSA-COSB)/(SINC*SINA) COSCS=(COSA*COSB-COSC)/(SINA*SINB) X(1,1)=AS*AS X(2,2)=BS*BS X(3,3)=CS*CS X(1,2)=2.*AS*BS*COSCS X(1,3)=2.*AS*CS*COSBS X(2,3)=2.*BS*CS*COSAS C CONVERT B(I,J) TO BETA(I,J) DO11I=1,N XL(I, 9)=XL(I, 9)*AS*AS*0.25 XL(I,10)=XL(I,10)*BS*BS*0.25 XL(I,11)=XL(I,11)*CS*CS*0.25 XL(I,12)=XL(I,12)*AS*BS*0.25 XL(I,13)=XL(I,13)*AS*CS*0.25 11 XL(I,14)=XL(I,14)*BS*CS*0.25 1 FORMAT(3F8.4,3F8.3) RETURN END SUBROUTINE RCELL(T,CC) C ******* C Subroutine to convert cell constants to real cell parameters. DIMENSION T(6),CC(6) RADS = 3.14159265359/180. AR = SQRT(T(1)) BR = SQRT(T(2)) CR = SQRT(T(3)) COSAR = T(4) / (2.0 * BR * CR) COSBR = T(5) / (2.0 * CR * AR) COSCR = T(6) / (2.0 * AR * BR) VR = AR * BR * CR * SQRT(1.0 - COSAR * COSAR - COSBR * COSBR - + COSCR * COSCR + 2.0 * COSAR * COSBR * COSCR) SINAR = SQRT(1.0 - COSAR * COSAR) SINBR = SQRT(1.0 - COSBR * COSBR) SINCR = SQRT(1.0 - COSCR * COSCR) A = BR * CR * SINAR / VR B = CR * AR * SINBR / VR C = AR * BR * SINCR / VR COSA =(COSBR * COSCR - COSAR) / (SINBR * SINCR) COSB =(COSCR * COSAR - COSBR) / (SINCR * SINAR) COSC =(COSAR * COSBR - COSCR) / (SINAR * SINBR) AL = ACOS(COSA) / RADS BE = ACOS(COSB) / RADS GA = ACOS(COSC) / RADS CC(1) = A CC(2) = B CC(3) = C CC(4) = AL CC(5) = BE CC(6) = GA RETURN C ******* END SUBROUTINE PROFIL C-----THIS SUBROUTINE ASSIGNS THE ADJUSTED VALUES TO THE PROFILE PARAM- C-----ETERS AND CALCULATES DIRECT CELL DIMENSIONS FROM THE CELL CONSTANT COMMON CUBIC,FORCED,HALPER,PAC,PUNCH,SKIP,STAR,TYPE,UNIAX, 1VERT,YES COMMON AH,B1,B2,BH,CH,DELTA,EPS,FACCOS,FACCUB,FACSIN,FLOG,FLOG2, 1FMAGDI,FMAGOB,FNUCDI,FNUCOB,OMEGA,PI,RELAXB,RELAXC,RELAXS,SCALE, 2SKEW,SOM,SUMCAL,SUMDFF,SUMDIF,SUMFOB,SUMOBS,SUMW,SUMWBS,SW,W,XN, 3YCALC,YOBS,ZERO,SLABDA,SLOG,ARG COMMON YALFA,ICENT,ICYCLE,YOMEG,IORD1,IORD2,IPONS,IPREV,IRL,YSTAP, 1KL,KM,KSTART,LCOORD,LCORD,LIM,LIMO,LLL,LMAT,LOUT,MAXCON,MAXS,MN,N, 2NCYCLE,NKEY,NLC,NN,NQC,NTEXT,NUM,YSS,NO COMMON ICODE(4889),ICON(4,9),JCON(4,9),LCONI(4,9),LCONJ(4,9),LN(4) 1,LP(71,14),MEQ(69),MTYP(69),NTYP(69),IOBS(10),ICALC(10) COMMON A(71,14),AL(3,3),ATEXT(69),B(6),C(4,98),CF(4),CONA(4,9), 1COSA(3),COSAF(3),COSAR(69,48),COSARI(48),DERIV(149),DERSTO(98,149) 2,DF(4),DUMP(98),DUMPA(500),F(2,21),FAC(4889),FESD(4889),FF(4889,2) 3,FLAG(71,14),FMAG(4889),FNUC(4889),FOBS(4889),H(48,3),HH(4889,3), 4HL(3),HN(3),HW(98),OCC(69),PAK(4889),PH(3),RM(48,9,9),RSA(69,3,3), 5RSB(69,3,3),S(2,21),SA(69),SB(69),SIDE(3),SINAR(69,48),SINARI(48), 6SJJ(4889),SM(48,3,3),SMM(149,149),SMOM(3),SNEX(69),SOME(98), 7SOMEGA(4889),T(48,3),TANN(98),TEMP(69),TEXT(20),THETA(4889),TR(48) 8,V(149),XL(71,14),SZ(71,14) COMMON NRD,NWD,NSD,NID,NFD,NDD,NCD,PHASEN(4889),SP(6,6),RFAC(6) 1,URABS,ABSOR LOGICAL HALPER,PAC,UNIAX,CUBIC,PUNCH,TYPE,STAR,SKIP, 1YES,VERT,FORCED DIMENSION AS(3),SINA(3) B1=XL(N+1,1) B2=XL(N+1,2) AH=XL(N+2,1) BH=XL(N+2,2) CH=XL(N+2,3) ZERO=XL(N+2,4) DO 1 I=1,3 1 AL(I,I)=XL(N+2,I+4) AL(2,3)=XL(N+2,8) AL(1,3)=XL(N+2,9) AL(1,2)=XL(N+2,10) T1=1. T2=1. T3=0. DO 2 I=1,3 J=1 IF(I.EQ.1)J=2 K=3 IF(I.EQ.3)K=2 AS(I)=SQRT(AL(I,I)) COSAF(I)=AL(J,K)/(2.*SQRT(AL(J,J)*AL(K,K))) SINA(I)=SQRT(1.-COSAF(I)*COSAF(I)) T1=T1*AL(I,I) T2=T2*AL(J,K) 2 T3=T3+AL(I,I)*AL(J,K)*AL(J,K) VL=SQRT(T1+0.25*(T2-T3)) DO 3 I=1,3 J=1 IF(I.EQ.1)J=2 K=3 IF(I.EQ.3)K=2 SIDE(I)=AS(J)*AS(K)*SINA(I)/VL 3 COSA(I)=(COSAF(J)*COSAF(K)-COSAF(I))/(SINA(J)*SINA(K)) RETURN END SUBROUTINE ITER C-----THIS SUBROUTINE IS CALLED ON EACH ITERATION COMMON CUBIC,FORCED,HALPER,PAC,PUNCH,SKIP,STAR,TYPE,UNIAX, 1VERT,YES COMMON AH,B1,B2,BH,CH,DELTA,EPS,FACCOS,FACCUB,FACSIN,FLOG,FLOG2, 1FMAGDI,FMAGOB,FNUCDI,FNUCOB,OMEGA,PI,RELAXB,RELAXC,RELAXS,SCALE, 2SKEW,SOM,SUMCAL,SUMDFF,SUMDIF,SUMFOB,SUMOBS,SUMW,SUMWBS,SW,W,XN, 3YCALC,YOBS,ZERO,SLABDA,SLOG,ARG COMMON YALFA,ICENT,ICYCLE,YOMEG,IORD1,IORD2,IPONS,IPREV,IRL,YSTAP, 1KL,KM,KSTART,LCOORD,LCORD,LIM,LIMO,LLL,LMAT,LOUT,MAXCON,MAXS,MN,N, 2NCYCLE,NKEY,NLC,NN,NQC,NTEXT,NUM,YSS,NO COMMON ICODE(4889),ICON(4,9),JCON(4,9),LCONI(4,9),LCONJ(4,9),LN(4) 1,LP(71,14),MEQ(69),MTYP(69),NTYP(69),IOBS(10),ICALC(10) COMMON A(71,14),AL(3,3),ATEXT(69),B(6),C(4,98),CF(4),CONA(4,9), 1COSA(3),COSAF(3),COSAR(69,48),COSARI(48),DERIV(149),DERSTO(98,149) 2,DF(4),DUMP(98),DUMPA(500),F(2,21),FAC(4889),FESD(4889),FF(4889,2) 3,FLAG(71,14),FMAG(4889),FNUC(4889),FOBS(4889),H(48,3),HH(4889,3), 4HL(3),HN(3),HW(98),OCC(69),PAK(4889),PH(3),RM(48,9,9),RSA(69,3,3), 5RSB(69,3,3),S(2,21),SA(69),SB(69),SIDE(3),SINAR(69,48),SINARI(48), 6SJJ(4889),SM(48,3,3),SMM(149,149),SMOM(3),SNEX(69),SOME(98), 7SOMEGA(4889),T(48,3),TANN(98),TEMP(69),TEXT(20),THETA(4889),TR(48) 8,V(149),XL(71,14),SZ(71,14) COMMON NRD,NWD,NSD,NID,NFD,NDD,NCD,PHASEN(4889),SP(6,6),RFAC(6) 1,URABS,ABSOR LOGICAL HALPER,PAC,UNIAX,CUBIC,PUNCH,TYPE,STAR,SKIP, 1YES,VERT,FORCED DIMENSION H1(3),H2(3),IN(3),I1(3),I2(3) SW=0. SUMW=0. SUMCAL=0. SUMOBS=0. SUMDIF=0. SUMWBS=0. SOM=0. SUMFOB=0. SUMDFF=0. FNUCDI=0. FNUCOB=0. FMAGDI=0. FMAGOB=0. NUM=0 K1=MAXS+MAXCON DO 1 I=1,K1 V(I)=0. DO 1 J=1,K1 1 SMM(I,J)=0. CALL CONPRE YES=.FALSE. IF(((IPONS.EQ.2).OR.(IPONS.EQ.3)).AND.(ICYCLE.EQ.NCYCLE))YES=.TRUE 1. TYPE=(LOUT.EQ.1).AND.(ICYCLE.EQ.NCYCLE) PUNCH=((IPONS.EQ.1).OR.(IPONS.EQ.3)).AND.(ICYCLE.EQ.NCYCLE) LLL=2 NK=NKEY NKEY=1 C DO 2 I=1,NK ! C 2 BACKSPACE NSD ! REWIND NSD ! A.HEWAT CHANGED FOR MAC 31/1/89 READ(NSD)DUMPA X=FLEES(0.) NN=FLEES(0.)+0.5 DO 3 I=1,3 3 X=FLEES(0.) DO 4 I=1,NN DO 5 J=1,5 5 X=FLEES(0.) SOMEGA(I)=0. FESD(I)=0. 4 FOBS(I)=0. CALL PROFIL YALFA=FLEES(0.) YSTAP=FLEES(0.) YOMEG=FLEES(0.) IF(.NOT.YES) GO TO 6 WRITE(NID,719) TEXT 719 FORMAT(20A4) TTSTEP=YSTAP/100. TTMIN=(YALFA-ZERO)/100. WRITE(NID,729) TTSTEP,SLABDA,TTMIN 729 FORMAT(6X,2H10,7X,1H0,F8.3,7X,1H1,7X,1H1,7X,1H0,F8.3/F8.3/) 6 IF(.NOT.TYPE)GOTO 13 WRITE(NWD,50)YALFA,YOMEG,YSTAP WRITE(NWD,9) 13 IPREV=1 NO=0 YSS=YALFA-YSTAP 1300 YSS=YSS+YSTAP Y=FLEES(0.) NO=NO+1 IF(ABS(Y).GE..001)GOTO 15 SOM=0. YOBS=0. YCALC=0. GOTO 108 15 YOBS=Y*SCALE NUM=NUM+1 W=FLEES(0.) SW=SW+W IORD1=FLEES(0.)+0.5 IORD2=FLEES(0.)+0.5 IF(IPREV.GT.IORD2)GOTO 40 IF(IPREV.LT.1) GO TO 40 !P.FISCHER 13 July 89 (Hewat) DO 21 J=IPREV,IORD2 JINDEX=J 21 CALL CALCUL(JINDEX) 40 IPREV=IORD2+1 108 CALL SUMMAT 8 IF(YSS.LE.(YOMEG-YSTAP)) GO TO 1300 IF(YES) WRITE(NID,739) 739 FORMAT(3X,5H-1000/2X,6H-10000) ARG=SQRT(3.14159265359)/(YSTAP*SCALE) DO 22 I=1,IORD2 X=SOMEGA(I) FOBS(I)=FOBS(I)/X IF(.NOT.((X.EQ.0.).OR.(FOBS(I).LT.0.)))GOTO 23 FOBS(I)=0. FESD(I)=0. 23 BB=SJJ(I)*ARG IF(TYPE)FESD(I)=SQRT(FESD(I)*250.)/X*BB FNUC(I)=FNUC(I)*BB FMAG(I)=FMAG(I)*BB 22 FOBS(I)=FOBS(I)*BB IF(.NOT.TYPE)GOTO 24 WRITE(NWD,25) K=IORD2/2 IF(KDIF.EQ.1) K=K+1 ! P Fischer, 20 Nov 86 ! DO 41 I=1,3 41 IN(I)=HN(I)+0.5 WRITE(NWD,27)((IN(I),I=1,3),J=1,2) DO 33 L=1,K I=L+K K1=FNUC(L)+FMAG(L)+0.5 HOLD=FOBS(L)-FLOAT(K1) L1=HOLD+SIGN(0.5,HOLD) M1=THETA(L)+ZERO+0.5 DO 29 J=1,3 HOLD=HH(L,J)*HN(J) 29 I1(J)=HOLD+SIGN(0.5,HOLD) IF(I-IORD2)30,30,32 30 K2=FNUC(I)+FMAG(I)+0.5 L2=FOBS(I)-FLOAT(K2)+0.5 M2=THETA(I)+ZERO+0.5 DO 31 J=1,3 HOLD=HH(I,J)*HN(J) 31 I2(J)=HOLD+SIGN(0.5,HOLD) GOTO 43 32 K2=0 L2=0 M2=0 FNUC(I)=0. FMAG(I)=0. FOBS(I)=0. FESD(I)=0. DO 34 J=1,3 34 I2(J)=0 43 LNUC=FNUC(L)+0.5 LMAG=FMAG(L)+0.5 LOBS=FOBS(L)+SIGN(0.5,FOBS(L)) LESD=FESD(L)+0.5 INUC=FNUC(I)+0.5 IMAG=FMAG(I)+0.5 KOBS=FOBS(I)+SIGN(0.5,FOBS(I)) IESD=FESD(I)+0.5 33 WRITE(NWD,35)(I1(J),J=1,3),M1,LNUC,LMAG,K1,LOBS,L1,LESD, 1(I2(J),J=1,3),M2,INUC,IMAG,K2,KOBS,L2,IESD 24 DO 36 I=1,IORD2 IF(FESD(I).GT.ABS(FNUC(I)+FMAG(I))) GO TO 3601 !A.W.HEWAT 28/9/89 SUMDFF=SUMDFF+ABS(FOBS(I)-ABS(FNUC(I)+FMAG(I))) SUMFOB=SUMFOB+FOBS(I) C...JKC DIV BY ZERO AJJ = FNUC(I)+FMAG(I) IF (AJJ.EQ.0.0) AJJ = 1.0E-10 X=FOBS(I)/AJJ FNUCDI=FNUCDI+ABS(FNUC(I)-FNUC(I)*X) FNUCOB=FNUCOB+FNUC(I)*X FMAGDI=FMAGDI+ABS(FMAG(I)-FMAG(I)*X) 3601 CONTINUE !A.W.HEWAT 28/9/89 36 FMAGOB=FMAGOB+FMAG(I)*X SKEW=0. 7 FORMAT(3I8/) 9 FORMAT(129H POSITION YOBS YCALC YOBS YCALC YOBS YCALC YOBS YCAHHHHHHHH 1LC YOBS YCALC YOBS YCALC YOBS YCALC YOBS YCALC YOBS YCALC YOHHHHHHHH 2BS YCALC) 14 FORMAT(1X,I6,5X) 17 FORMAT(/) 25 FORMAT(129H0 H K L POS INUC IMAG ITOT IOBS ITER1630 1DIF ESD * H K L POS INUC IMAG ITOT IOBS ITER1640 2DIF ESD) ITER1650 27 FORMAT(1X,I3,2I4,56X,I2,2I3) ITER1660 35 FORMAT(1X,I3,2I4,I7,4I8,I7,I5,4H * ,I3,2I4,I7,4I8,I7,I5) ITER1670 50 FORMAT(16H0DIAGRAM IS FROM,F8.1,1X,2HTO,F8.1,1X,11HIN STEPS OF, AF5.1,1X,16H*1/100THS DEGREE) ITER1690 RETURN ITER1700 END ITER1710 SUBROUTINE CONPRE ITER1620 C-----THIS SUBROUTINE PREPARES THE ELEMENTS OF THE BORDERING MATRIX DUE C-----TO THE CONSTRAINT FUNCTIONS. COMMON CUBIC,FORCED,HALPER,PAC,PUNCH,SKIP,STAR,TYPE,UNIAX, 1VERT,YES COMMON AH,B1,B2,BH,CH,DELTA,EPS,FACCOS,FACCUB,FACSIN,FLOG,FLOG2, 1FMAGDI,FMAGOB,FNUCDI,FNUCOB,OMEGA,PI,RELAXB,RELAXC,RELAXS,SCALE, 2SKEW,SOM,SUMCAL,SUMDFF,SUMDIF,SUMFOB,SUMOBS,SUMW,SUMWBS,SW,W,XN, 3YCALC,YOBS,ZERO,SLABDA,SLOG,ARG COMMON YALFA,ICENT,ICYCLE,YOMEG,IORD1,IORD2,IPONS,IPREV,IRL,YSTAP, 1KL,KM,KSTART,LCOORD,LCORD,LIM,LIMO,LLL,LMAT,LOUT,MAXCON,MAXS,MN,N, 2NCYCLE,NKEY,NLC,NN,NQC,NTEXT,NUM,YSS,NO COMMON ICODE(4889),ICON(4,9),JCON(4,9),LCONI(4,9),LCONJ(4,9),LN(4) 1,LP(71,14),MEQ(69),MTYP(69),NTYP(69),IOBS(10),ICALC(10) COMMON A(71,14),AL(3,3),ATEXT(69),B(6),C(4,98),CF(4),CONA(4,9), 1COSA(3),COSAF(3),COSAR(69,48),COSARI(48),DERIV(149),DERSTO(98,149) 2,DF(4),DUMP(98),DUMPA(500),F(2,21),FAC(4889),FESD(4889),FF(4889,2) 3,FLAG(71,14),FMAG(4889),FNUC(4889),FOBS(4889),H(48,3),HH(4889,3), 4HL(3),HN(3),HW(98),OCC(69),PAK(4889),PH(3),RM(48,9,9),RSA(69,3,3), 5RSB(69,3,3),S(2,21),SA(69),SB(69),SIDE(3),SINAR(69,48),SINARI(48), 6SJJ(4889),SM(48,3,3),SMM(149,149),SMOM(3),SNEX(69),SOME(98), 7SOMEGA(4889),T(48,3),TANN(98),TEMP(69),TEXT(20),THETA(4889),TR(48) 8,V(149),XL(71,14),SZ(71,14) COMMON NRD,NWD,NSD,NID,NFD,NDD,NCD,PHASEN(4889),SP(6,6),RFAC(6) 1,URABS,ABSOR LOGICAL HALPER,PAC,UNIAX,CUBIC,PUNCH,TYPE,STAR,SKIP, 1YES,VERT,FORCED IF(MAXCON.EQ.0)GOTO 1 DO 2 I=1,MAXCON DO 2 J=1,MAXS 2 C(I,J)=0. DO 3 I=1,MAXCON CF(I)=DF(I) LNI=LN(I) DO 3 J=1,LNI IL=ICON(I,J) JL=JCON(I,J) AA=CONA(I,J) IF(I.GT.NLC)GOTO 4 LPIJ=LP(IL,JL) C(I,LPIJ)=AA*A(IL,JL)+C(I,LPIJ) CF(I)=CF(I)-AA*XL(IL,JL) GOTO 3 4 X=XL(IL,JL) II=LCONI(I,J) JJ=LCONJ(I,J) Y=XL(II,JJ) LPIJ=LP(IL,JL) C(I,LPIJ)=AA*A(IL,JL)*Y+C(I,LPIJ) LPIJ=LP(II,JJ) C(I,LPIJ)=AA*A(II,JJ)*X+C(I,LPIJ) CF(I)=CF(I)-AA*X*Y*0.5 3 CONTINUE 1 RETURN END SUBROUTINE SUMMAT C-----THIS SUBROUTINE FORMS THE ELEMENTS OF THE NORMAL MATRIX AND VECTOR C-----AND CALCULATES SEVERAL SUMS COMMON CUBIC,FORCED,HALPER,PAC,PUNCH,SKIP,STAR,TYPE,UNIAX, 1VERT,YES COMMON AH,B1,B2,BH,CH,DELTA,EPS,FACCOS,FACCUB,FACSIN,FLOG,FLOG2, 1FMAGDI,FMAGOB,FNUCDI,FNUCOB,OMEGA,PI,RELAXB,RELAXC,RELAXS,SCALE, 2SKEW,SOM,SUMCAL,SUMDFF,SUMDIF,SUMFOB,SUMOBS,SUMW,SUMWBS,SW,W,XN, 3YCALC,YOBS,ZERO,SLABDA,SLOG,ARG COMMON YALFA,ICENT,ICYCLE,YOMEG,IORD1,IORD2,IPONS,IPREV,IRL,YSTAP, 1KL,KM,KSTART,LCOORD,LCORD,LIM,LIMO,LLL,LMAT,LOUT,MAXCON,MAXS,MN,N, 2NCYCLE,NKEY,NLC,NN,NQC,NTEXT,NUM,YSS,NO COMMON ICODE(4889),ICON(4,9),JCON(4,9),LCONI(4,9),LCONJ(4,9),LN(4) 1,LP(71,14),MEQ(69),MTYP(69),NTYP(69),IOBS(10),ICALC(10) COMMON A(71,14),AL(3,3),ATEXT(69),B(6),C(4,98),CF(4),CONA(4,9), 1COSA(3),COSAF(3),COSAR(69,48),COSARI(48),DERIV(149),DERSTO(98,149) 2,DF(4),DUMP(98),DUMPA(500),F(2,21),FAC(4889),FESD(4889),FF(4889,2) 3,FLAG(71,14),FMAG(4889),FNUC(4889),FOBS(4889),H(48,3),HH(4889,3), 4HL(3),HN(3),HW(98),OCC(69),PAK(4889),PH(3),RM(48,9,9),RSA(69,3,3), 5RSB(69,3,3),S(2,21),SA(69),SB(69),SIDE(3),SINAR(69,48),SINARI(48), 6SJJ(4889),SM(48,3,3),SMM(149,149),SMOM(3),SNEX(69),SOME(98), 7SOMEGA(4889),T(48,3),TANN(98),TEMP(69),TEXT(20),THETA(4889),TR(48) 8,V(149),XL(71,14),SZ(71,14) COMMON NRD,NWD,NSD,NID,NFD,NDD,NCD,PHASEN(4889),SP(6,6),RFAC(6) 1,URABS,ABSOR LOGICAL HALPER,PAC,UNIAX,CUBIC,PUNCH,TYPE,STAR,SKIP, SUMM 230 1YES,VERT,FORCED SUMM 240 LOGICAL VERT1 CJH DIMENSION LPM(10) SUMM 250 IF(YOBS.EQ.0.)GOTO 20 SUMM 260 DO 1 J=1,10 SUMM 270 1 LPM(J)=LP(N+2,J) SUMM 280 LP2=LP(N+1,2) SUMM 290 LP4=LP(N+2,4) SUMM 300 DO 2 J=1,MAXS SUMM 310 2 DERIV(J)=0. SUMM 320 VERT1=LIM.LT.0 CJH YCALC=0. SUMM 330 IL=0 SUMM 340 C-----CALCULATE THE CONTRIBUTION OF THE REFLECTIONS ORD1 TO ORD2 TO THE SUMM 350 C-----DERIVATIVES W.R.T. THE PROFILE INTENSITY YOBS SUMM 360 DO 3 J=IORD1,IORD2 SUMM 370 IL=IL+1 SUMM 380 J10=MOD(J,LIMO)+1 SUMM 390 VERT=THETA(J).LE.FLOAT(LIM) SUMM 400 DELTA=YSS-THETA(J)-ZERO SUMM 410 DELT=DELTA*DELTA SUMM 420 TL=HW(J10) SUMM 430 BB=TL*TL SUMM 440 IF(.NOT.VERT)GOTO 4 SUMM 450 Y=DELT*SIGN(1.,DELTA) SUMM 460 Z=1.-B2*Y/TANN(J10) SUMM 470 IF(Z.LE.0.)Z=.0001 SUMM 480 GOTO 5 SUMM 490 4 Z=1. SUMM 500 IF(.NOT.VERT1) GO TO 5 CJH COT2TH=0.5*(1.-TANN(J10)*TANN(J10))/TANN(J10) CJH ASYM=B2*COT2TH CJH 5 PHI=EXP(-DELT/BB)/TL CJH IF(.NOT.VERT1) GO TO 50 CJH EXP1=EXP(-(DELTA/TL)**2) CJH EXP2=EXP(-((DELTA+ASYM/16.0)/TL)**2) CJH EXP3=EXP(-((DELTA+ASYM/4.0)/TL)**2) CJH EXP4=EXP(-((DELTA+9.0*ASYM/16.0)/TL)**2) CJH EXP5=EXP(-((DELTA+ASYM)/TL)**2) CJH PHI=(EXP1+4.0*EXP2+2.0*EXP3+4.0*EXP4+EXP5)/(12.0*TL) CJH DERHW=(DELTA**2*EXP1 + 4.0*(DELTA+ASYM/16.0)**2*EXP2 CJH 1 + 2.0*(DELTA+ASYM/4.)**2*EXP3 + 4.0*(DELTA+9.*ASYM/16.)**2*EXP4 CJH 1 + (DELTA+ASYM)**2*EXP5)/(6.0*TL**3*PHI) - 1.0 CJH DERZ=(DELTA*EXP1 + 4.0*(DELTA+ASYM/16.0)*EXP2 CJH 1 + 2.0*(DELTA+ASYM/4.0)*EXP3 + 4.0*(DELTA+9.*ASYM/16.0)*EXP4 CJH 2 + (DELTA+ASYM)*EXP5) / (12.0*TL*PHI) CJH DERASY=((DELTA+ASYM/16.0)*EXP2 + 2.0*(DELTA+ASYM/4.0)*EXP3 CJH 1 + 9.0*(DELTA+9.*ASYM/16.0)*EXP4 + 4.0*(DELTA+ASYM)*EXP5) CJH 2 / (24.0*TL**3*PHI) CJH 50 OMEGA=SJJ(J)*PHI*Z CJH FNN=FNUC(J)+FMAG(J) SUMM 520 YCALC=YCALC+OMEGA*FNN SUMM 530 DUMP(IL)=OMEGA*FNN SUMM 540 SOME(IL)=OMEGA SUMM 550 X=2.*DELT/BB-1. SUMM 560 DO 3 K=1,MAXS SUMM 570 DER=1. SUMM 580 DO 6 M=1,3 SUMM 590 IF(LPM(M).EQ.K)DER=X SUMM 600 IF((LPM(M).EQ.K).AND.VERT1)DER=DERHW CJH 6 CONTINUE SUMM 610 X1=0. SUMM 620 C I THINK THE FOLLOWING LACKS A FACTOR 1/TANN(J10) CJH IF(VERT)X1=B2*SIGN(1.,DELTA)*BB/Z SUMM 630 IF(LP4.EQ.K)DER=DELTA*(1.+X1) SUMM 640 IF((LP4.EQ.K).AND.VERT1) DER=DERZ CJH IF((LP2.EQ.K).AND.VERT)DER=Y/Z SUMM 650 IF((LP2.EQ.K).AND.VERT1) DER=DERASY CJH DO 7 M=5,10 SUMM 660 C I THINK THE FOLLOWING LACKS A FACTOR (1.+X1) CJH IF(LPM(M).EQ.K)DER=DELTA SUMM 670 IF((LPM(M).EQ.K).AND.VERT1) DER=DERZ CJHTEMP C IF((LPM(M).EQ.K).AND.VERT1)DER=DERZ+0.5*BB*B2*(1.+COT2TH**2) CJH C 1 *DERASY CJH 7 CONTINUE SUMM 680 3 DERIV(K)=DERSTO(J10,K)*DER*OMEGA+DERIV(K) SUMM 690 C-----FORM SUMS SUMM 700 YCALC=YCALC/SCALE SUMM 710 YOBS=YOBS/SCALE SUMM 720 DELTA=YOBS-YCALC SUMM 730 DO 8 J=IORD1,IORD2 SUMM 740 K=J-IORD1+1 SUMM 750 C Divide by zero! IF (YCALC.EQ.0.0) YCALC = 1.0E-10 X=DUMP(K)/YCALC SUMM 760 IF(TYPE)FESD(J)=X*X/W+FESD(J) SUMM 770 FOBS(J)=FOBS(J)+X*YOBS SUMM 780 8 SOMEGA(J)=SOMEGA(J)+SOME(K) SUMM 790 SUMW=SUMW+W*DELTA*DELTA SUMM 800 SUMWBS=SUMWBS+W*YOBS*YOBS SUMM 810 SUMCAL=SUMCAL+YCALC SUMM 820 SUMOBS=SUMOBS+ABS(YOBS) SUMM 830 SUMDIF=SUMDIF+ABS(DELTA) SUMM 840 DELTA=DELTA*SCALE SUMM 850 C-----SUM NORMAL MATRIX ELEMENTS AND VECTOR SUMM 860 DO 9 J=1,MAXS SUMM 870 X=W*DERIV(J) SUMM 880 V(J)=V(J)+DELTA*X SUMM 890 DO 9 K=J,MAXS SUMM 900 9 SMM(J,K)=SMM(J,K)+X*DERIV(K) SUMM 910 20 IOBS(NO)=YOBS+SIGN(0.5,YOBS) SUMM 920 ICALC(NO)=YCALC+0.5 SUMM 930 IF(.NOT.((NO.EQ.10).OR.(YSS.GE.YOMEG)))GOTO 10 SUMM 940 K=IFIX(YSS-FLOAT(NO-1)*YSTAP) SUMM 950 C-----PRINT CALCULATED AND OBSERVED PROFILE INTENSITIES SUMM 960 IF(TYPE)WRITE(NWD,12)K,(IOBS(J),ICALC(J),J=1,NO) SUMM 970 C-----PUNCH CALCULATED AND OBSERVED PROFILE INTENSITIES SUMM 980 IF(YES)WRITE(NID,14)(IOBS(J),J=1,NO) SUMM 990 IF(YES)WRITE(NID,14)(ICALC(J),J=1,NO) NO=0 SUMM1000 10 CONTINUE SUMM1010 12 FORMAT(1X,I5,1X,20I6) SUMM1020 14 FORMAT(10I8) SUMM1030 RETURN SUMM1040 END SUMM1050 SUBROUTINE MATRIX MATR 0 C-----THIS SUBROUTING BORDERS THE NORMAL MATRIX WITH THE CONSTRAINT FUNCMATR 10 C-----TIONS,INVERTS THIS ENLARGED MATRIX AND SOLVES THE SYSTEM MATR 20 COMMON CUBIC,FORCED,HALPER,PAC,PUNCH,SKIP,STAR,TYPE,UNIAX, MATR 30 1VERT,YES MATR 40 COMMON AH,B1,B2,BH,CH,DELTA,EPS,FACCOS,FACCUB,FACSIN,FLOG,FLOG2, MATR 50 1FMAGDI,FMAGOB,FNUCDI,FNUCOB,OMEGA,PI,RELAXB,RELAXC,RELAXS,SCALE, MATR 60 2SKEW,SOM,SUMCAL,SUMDFF,SUMDIF,SUMFOB,SUMOBS,SUMW,SUMWBS,SW,W,XN, MATR 70 3YCALC,YOBS,ZERO,SLABDA,SLOG,ARG MATR 80 COMMON YALFA,ICENT,ICYCLE,YOMEG,IORD1,IORD2,IPONS,IPREV,IRL,YSTAP,MATR 90 1KL,KM,KSTART,LCOORD,LCORD,LIM,LIMO,LLL,LMAT,LOUT,MAXCON,MAXS,MN,N,MATR 100 2NCYCLE,NKEY,NLC,NN,NQC,NTEXT,NUM,YSS,NO MATR 110 COMMON ICODE(4889),ICON(4,9),JCON(4,9),LCONI(4,9),LCONJ(4,9),LN(4) 1,LP(71,14),MEQ(69),MTYP(69),NTYP(69),IOBS(10),ICALC(10) COMMON A(71,14),AL(3,3),ATEXT(69),B(6),C(4,98),CF(4),CONA(4,9), 1COSA(3),COSAF(3),COSAR(69,48),COSARI(48),DERIV(149),DERSTO(98,149) 2,DF(4),DUMP(98),DUMPA(500),F(2,21),FAC(4889),FESD(4889),FF(4889,2) 3,FLAG(71,14),FMAG(4889),FNUC(4889),FOBS(4889),H(48,3),HH(4889,3), 4HL(3),HN(3),HW(98),OCC(69),PAK(4889),PH(3),RM(48,9,9),RSA(69,3,3), 5RSB(69,3,3),S(2,21),SA(69),SB(69),SIDE(3),SINAR(69,48),SINARI(48), 6SJJ(4889),SM(48,3,3),SMM(149,149),SMOM(3),SNEX(69),SOME(98), 7SOMEGA(4889),T(48,3),TANN(98),TEMP(69),TEXT(20),THETA(4889),TR(48) 8,V(149),XL(71,14),SZ(71,14) COMMON NRD,NWD,NSD,NID,NFD,NDD,NCD,PHASEN(4889),SP(6,6),RFAC(6) 1,URABS,ABSOR LOGICAL HALPER,PAC,UNIAX,CUBIC,PUNCH,TYPE,STAR,SKIP, MATR 230 1YES,VERT,FORCED MATR 240 DIMENSION KS(20) MATR 250 DO 1 I=1,MAXS MATR 260 DO 1 J=I,MAXS MATR 270 1 SMM(J,I)=SMM(I,J) MATR 280 C-----BORDER THE MATRIX AND ENLARGE THE VECTOR MATR 290 IF(MAXCON)13,13,14 MATR 300 14 DO 2 I=1,MAXCON MATR 310 IMAXS=I+MAXS MATR 320 V( IMAXS)=CF(I)*SW MATR 330 DO 2 J=1,MAXS MATR 340 SMM(J, IMAXS)=C(I,J)*SW MATR 350 2 SMM( IMAXS,J)=SMM(J, IMAXS) MATR 360 13 SKEW=DTSLIN(MAXS+MAXCON) MATR 370 OOM=SUMW/FLOAT(NUM-MAXS+MAXCON)*(SCALE*SCALE) MATR 380 C-----CALCULATE THE COVARIANCE MATRIX MATR 390 DO 3 I=1,MAXS MATR 400 DO 3 J=1,MAXS MATR 410 3 SMM(I,J)=SMM(I,J)*OOM MATR 420 C-----TEST ON CONVERGENCE MATR 430 DO 4 I=1,MAXS MATR 440 IF(SMM(I,I).LT.0.) WRITE(*,499) I,I,SMM(I,I) 499 FORMAT(' ATTENTION. SMM(',I3,',',I3,') NEGATIVE=',F10.4) X1=SQRT(ABS(SMM(I,I)))*EPS ! A.H. 18 JAN'88 MATR 450 IF(ABS(V(I)).GT.X1)GOTO 5 MATR 460 4 CONTINUE MATR 470 IF(ICYCLE.EQ.NCYCLE)GOTO 5 MATR 480 ICYCLE=NCYCLE-1 MATR 490 FORCED=.TRUE. MATR 500 C-----PRINT CORRELATION MATRIX MATR 510 5 IF(.NOT.((LMAT.EQ.1).AND.(ICYCLE.EQ.NCYCLE)))GOTO 6 MATR 520 WRITE(NWD,7) MATR 530 IA=1 MATR 540 LIM=MAXS MATR 550 15 IF((LIM-IA+1).GT.20) LIM=IA+19 MATR 560 WRITE(NWD,8)(I,I=IA,LIM) MATR 570 DO 11 I=1,MAXS MATR 580 L=0 MATR 590 X=SMM(I,I) MATR 600 DO 16 J=IA,LIM MATR 610 L=L+1 MATR 620 Y=SMM(J,J)*X MATR 630 16 KS(L)=100.*SMM(I,J)/(SQRT(ABS(Y))*SIGN(1.,Y))+0.5 MATR 640 11 WRITE(NWD,12)I,(KS(J),J=1,L) MATR 650 IF(LIM.GE.MAXS)GOTO 6 MATR 660 IA=LIM+1 MATR 670 LIM=MAXS MATR 680 GOTO 15 MATR 690 6 CONTINUE MATR 700 7 FORMAT(20H0CORRELATION MATRIX=) MATR 710 8 FORMAT(5H0 ,20I6) MATR 720 12 FORMAT(I5,20I6) MATR 730 RETURN MATR 740 END MATR 750 FUNCTION DTSLIN(N) DTSL 0 C C-----DTSLIN SUBROUTINE SOLVES THE MATRIX EQUATION A.X=B FOR X AND STOREDTSL 10 C-----THE RESULT X IN B.THE SIZE OF THE SYSTEM IS N,WHERE N.LE.149 DTSL 20 DIMENSION V(149),M(149),U(149) DTSL 30 COMMON CUBIC,FORCED,HALPER,PAC,PUNCH,SKIP,STAR,TYPE,UNIAX, DTSL 40 1VERT,YES DTSL 50 COMMON AH,B1,B2,BH,CH,DELTA,EPS,FACCOS,FACCUB,FACSIN,FLOG,FLOG2, DTSL 60 1FMAGDI,FMAGOB,FNUCDI,FNUCOB,OMEGA,PI,RELAXB,RELAXC,RELAXS,SCALE, DTSL 70 2SKEW,SOM,SUMCAL,SUMDFF,SUMDIF,SUMFOB,SUMOBS,SUMW,SUMWBS,SW,W,XN, DTSL 80 3YCALC,YOBS,ZERO,SLABDA,SLOG,ARG DTSL 90 COMMON YALFA,ICENT,ICYCLE,YOMEG,IORD1,IORD2,IPONS,IPREV,IRL,YSTAP,DTSL 100 1KL,KM,KSTART,LCOORD,LCORD,LIM,LIMO,LLL,LMAT,LOUT,MAXCON,MAXS,MN,L,DTSL 110 2NCYCLE,NKEY,NLC,NN,NQC,NTEXT,NUM,YSS,NO DTSL 120 COMMON ICODE(4889),ICON(4,9),JCON(4,9),LCONI(4,9),LCONJ(4,9),LN(4) 1,LP(71,14),MEQ(69),MTYP(69),NTYP(69),IOBS(10),ICALC(10) COMMON Q(71,14),AL(3,3),ATEXT(69),G(6),Z(4,98),CF(4),CONA(4,9), 1COSA(3),COSAF(3),COSAR(69,48),COSARI(48),DERIV(149),DERSTO(98,149) 2,DF(4),DUMP(98),DUMPA(500),F(2,21),FAC(4889),FESD(4889),FF(4889,2) 3,FLAG(71,14),FMAG(4889),FNUC(4889),FOBS(4889),H(48,3),HH(4889,3), 4HL(3),HN(3),HW(98),OCC(69),PAK(4889),PH(3),RM(48,9,9),RSA(69,3,3), 5RSB(69,3,3),X(2,21),SA(69),SB(69),SIDE(3),SINAR(69,48),SINARI(48), 6SJJ(4889),SM(48,3,3),A (149,149),SMOM(3),SNEX(69),SOME(98), 7SOMEGA(4889),T(48,3),TANN(98),TEMP(69),TEXT(20),THETA(4889),TR(48) 8,B(149),XL(71,14),SZ(71,14) COMMON NRD,NWD,NSD,NID,NFD,NDD,NCD,PHASEN(4889),SP(6,6),RFAC(6) 1,URABS,ABSOR LOGICAL HALPER,PAC,UNIAX,CUBIC,PUNCH,TYPE,STAR,SKIP, DTSL 240 1YES,VERT,FORCED DTSL 250 I=1 DTSL 260 3 C=0. DTSL 270 D=0. DTSL 280 J=1 DTSL 290 2 IF(J.GT.N)GOTO 1 DTSL 300 C=A(I,J)*A(I,J)+C DTSL 310 D=A(J,I)*A(J,I)+D DTSL 320 J=J+1 DTSL 330 GOTO 2 DTSL 340 1 V(I)=SQRT(C) DTSL 350 U(I)=SQRT(D) DTSL 360 I=I+1 DTSL 370 IF(I.LE.N)GOTO 3 DTSL 380 I=1 DTSL 390 41 C=1./U(I) DTSL 400 J=1 DTSL 410 40 A(J,I)=A(J,I)*C DTSL 420 J=J+1 DTSL 430 IF(J.LE.N)GOTO 40 DTSL 440 I=I+1 DTSL 450 IF(I.LE.N)GOTO 41 DTSL 460 E=1. DTSL 470 K=1 DTSL 480 15 R=-1. DTSL 490 I=K DTSL 500 8 C=0. DTSL 510 J=1 DTSL 520 K1=K-1 DTSL 530 5 IF(J.GT.K1)GOTO 4 DTSL 540 C=A(I,J)*A(J,K)+C DTSL 550 J=J+1 DTSL 560 GOTO 5 DTSL 570 4 A(I,K)=A(I,K)-C DTSL 580 S=ABS(A(I,K)/V(I)) DTSL 590 IF(S-R)6,6,7 DTSL 600 7 R=S DTSL 610 M(K)=I DTSL 620 6 I=I+1 DTSL 630 IF(I.LE.N)GOTO 8 DTSL 640 KP=M(K) DTSL 650 V(KP)=V(K) DTSL 660 J=1 DTSL 670 14 R=A(K,J) DTSL 680 IF(J-K)9,9,10 DTSL 690 9 KP=M(K) DTSL 700 A(K,J)=A(KP,J) DTSL 710 GOTO 11 DTSL 720 10 C=0. DTSL 730 I=1 DTSL 740 K1=K-1 DTSL 750 13 IF(I.GT.K1)GOTO 12 DTSL 760 C=A(K,I)*A(I,J)+C DTSL 770 I=I+1 DTSL 780 GOTO 13 DTSL 790 12 KP=M(K) DTSL 800 A(K,J)=(A(KP,J)-C)/A(K,K) DTSL 810 11 KP=M(K) DTSL 820 IF(KP.NE.K)A(KP,J)=-R DTSL 830 J=J+1 DTSL 840 IF(J.LE.N)GOTO 14 DTSL 850 E=A(K,K)*E DTSL 860 K=K+1 DTSL 870 IF(K.LE.N)GOTO 15 DTSL 880 DTSLIN=E DTSL 890 K=1 DTSL 900 18 R=B(K) DTSL 910 C=0. DTSL 920 I=1 DTSL 930 K1=K-1 DTSL 940 17 IF(I.GT.K1)GOTO 16 DTSL 950 C=A(K,I)*B(I)+C DTSL 960 I=I+1 DTSL 970 GOTO 17 DTSL 980 16 KP=M(K) DTSL 990 B(K)=(B(KP)-C)/A(K,K) DTSL1000 IF(KP.NE.K)B(KP)=-R DTSL1010 K=K+1 DTSL1020 IF(K.LE.N)GOTO 18 DTSL1030 K=N DTSL1040 21 C=0. DTSL1050 I=K+1 DTSL1060 20 IF(I.GT.N)GOTO 19 DTSL1070 C=A(K,I)*B(I)+C DTSL1080 I=I+1 DTSL1090 GOTO 20 DTSL1100 19 B(K)=B(K)-C DTSL1110 K=K-1 DTSL1120 IF(K.GE.1)GOTO 21 DTSL1130 K=N DTSL1140 31 J=K+1 DTSL1150 23 IF(J.GT.N)GOTO 22 DTSL1160 V(J)=A(K,J) DTSL1170 A(K,J)=0. DTSL1180 J=J+1 DTSL1190 GOTO 23 DTSL1200 22 A(K,K)=1./A(K,K) DTSL1210 J=K-1 DTSL1220 27 IF(J.LT.1)GOTO 24 DTSL1230 C=0. DTSL1240 I=J+1 DTSL1250 26 IF(I.GT.K)GOTO 25 DTSL1260 C=A(K,I)*A(I,J)+C DTSL1270 I=I+1 DTSL1280 GOTO 26 DTSL1290 25 A(K,J)=-C/A(J,J) DTSL1300 J=J-1 DTSL1310 IF(J.GE.1)GOTO 27 DTSL1320 24 J=1 DTSL1330 30 C=0. DTSL1340 I=K+1 DTSL1350 29 IF(I.GT.N)GOTO 28 DTSL1360 C=V(I)*A(I,J)+C DTSL1370 I=I+1 DTSL1380 GOTO 29 DTSL1390 28 A(K,J)=A(K,J)-C DTSL1400 J=J+1 DTSL1410 IF(J.LE.N)GOTO 30 DTSL1420 K=K-1 DTSL1430 IF(K.GE.1)GOTO 31 DTSL1440 K=N DTSL1450 34 IF(M(K).EQ.K)GOTO 32 DTSL1460 I=1 DTSL1470 33 R=A(I,K) DTSL1480 KP=M(K) DTSL1490 A(I,K)=-A(I,KP) DTSL1500 A(I,KP)=R DTSL1510 I=I+1 DTSL1520 IF(I.LE.N)GOTO 33 DTSL1530 32 K=K-1 DTSL1540 IF(K.GE.1)GOTO 34 DTSL1550 I=1 DTSL1560 43 B(I)=B(I)/U(I) DTSL1570 J=1 DTSL1580 42 A(I,J)=A(I,J)/U(I) DTSL1590 J=J+1 DTSL1600 IF(J.LE.N)GOTO 42 DTSL1610 I=I+1 DTSL1620 IF(I.LE.N)GOTO 43 DTSL1630 RETURN DTSL1640 END DTSL1650 SUBROUTINE DIRECT DIRE 0 C-----THIS SUBROUTINE PRINTS AND CALCULATES THE DIRECT CELL DIMENSIONS DIRE 10 C-----AND THEIR E.S.D. FROM THE CELL CONSTANTS AND THEIR E.S.D. DIRE 20 COMMON CUBIC,FORCED,HALPER,PAC,PUNCH,SKIP,STAR,TYPE,UNIAX, DIRE 30 1VERT,YES DIRE 40 COMMON AH,B1,B2,BH,CH,DELTA,EPS,FACCOS,FACCUB,FACSIN,FLOG,FLOG2, DIRE 50 1FMAGDI,FMAGOB,FNUCDI,FNUCOB,OMEGA,PI,RELAXB,RELAXC,RELAXS,SCALE, DIRE 60 2SKEW,SOM,SUMCAL,SUMDFF,SUMDIF,SUMFOB,SUMOBS,SUMW,SUMWBS,SW,W,XN, DIRE 70 3YCALC,YOBS,ZERO,SLABDA,SLOG,ARG DIRE 80 COMMON YALFA,ICENT,ICYCLE,YOMEG,IORD1,IORD2,IPONS,IPREV,IRL,YSTAP,DIRE 90 1KL,KM,KSTART,LCOORD,LCORD,LIM,LIMO,LLL,LMAT,LOUT,MAXCON,MAXS,MN,N,DIRE 100 2NCYCLE,NKEY,NLC,NN,NQC,NTEXT,NUM,YSS,NO DIRE 110 COMMON ICODE(4889),ICON(4,9),JCON(4,9),LCONI(4,9),LCONJ(4,9),LN(4) 1,LP(71,14),MEQ(69),MTYP(69),NTYP(69),IOBS(10),ICALC(10) COMMON A(71,14),AL(3,3),ATEXT(69),B(6),C(4,98),CF(4),CONA(4,9), 1COSA(3),COSAF(3),COSAR(69,48),COSARI(48),DEROV(149),DERSTO(98,149) 2,DF(4),DUMP(98),DUMPA(500),F(2,21),FAC(4889),FESD(4889),FF(4889,2) 3,FLAG(71,14),FMAG(4889),FNUC(4889),FOBS(4889),H(48,3),HH(4889,3), 4HL(3),HN(3),HW(98),OCC(69),PAK(4889),PH(3),RM(48,9,9),RSA(69,3,3), 5RSB(69,3,3),R(2,21),SO(69),SB(69),SIDE(3),SINAR(69,48),SINARI(48), 6SJJ(4889),SQ(48,3,3),SMM(149,149),SMOM(3),SNEX(69),SOME(98), 7SOMEGA(4889),T(48,3),TANN(98),TEMP(69),TEXT(20),THETA(4889),TR(48) 8,Q(149),XL(71,14),SZ(71,14) COMMON NRD,NWD,NSD,NID,NFD,NDD,NCD,PHASEN(4889),SM(6,6),RFAC(6) 1,URABS,ABSOR LOGICAL HALPER,PAC,UNIAX,CUBIC,PUNCH,TYPE,STAR,SKIP, DIRE 230 1YES,VERT,FORCED DIRE 240 DIMENSION V(6), DERIV(6),SA(6),SSA(6) DIRE 250 DIMENSION BIJ(6),SBIJ(6),FAK(6) RAD=3.14159265359/360. DIRE 260 DO 1 J=1,6 DIRE 270 V(J)=XL(N+2,J+4) DIRE 280 K=LP(N+2,J+4) DIRE 290 X=A(N+2,J+4) DIRE 300 IF(K.EQ.0)X=0. DIRE 310 DO 1 L=J,6 DIRE 320 M=LP(N+2,L+4) DIRE 330 SM(L,J)=0. DIRE 340 IF((M.EQ.0).OR.(K.EQ.0))GOTO 1 DIRE 350 SM(L,J)=SMM(M,K)*X*A(N+2,L+4) DIRE 360 1 SM(J,L)=SM(L,J) DIRE 370 VOLUME=V(1)*V(2)*V(3)+.25*V(4)*V(5)*V(6) DIRE 380 DO 3 I=1,3 DIRE 390 VOLUME=VOLUME-V(I)*V(I+3)*V(I+3)*.25 DIRE 400 J=I+1 DIRE 410 IF(J.GT.3)J=J-3 DIRE 420 K=J+1 DIRE 430 IF(K.GT.3)K=K-3 DIRE 440 DERIV(I)=V(J)*V(K)-V(I+3)**2*.25 DIRE 450 3 DERIV(I+3)=.25*(V(J+3)*V(K+3)-V(I)*V(I+3)*2.) DIRE 460 SVOL=0. DIRE 470 DO 4 I=1,6 DIRE 480 DO 4 J=I,6 DIRE 490 X=1. DIRE 500 IF(I.NE.J)X=2. DIRE 510 4 SVOL=SVOL+SM(I,J)*DERIV(I)*DERIV(J)*X DIRE 520 DO 5 I=1,3 DIRE 530 DO 6 J=1,6 DIRE 540 6 DERIV(J)=0. DIRE 550 J=I+1 DIRE 560 IF(J.GT.3)J=J-3 DIRE 570 K=J+1 DIRE 580 IF(K.GT.3)K=K-3 DIRE 590 SA(I)=SQRT(V(I)) DIRE 600 SA(I+3)=V(I+3)/SQRT(V(J)*V(K))*0.5 DIRE 610 SSA(I)=SM(I,I)/(4.*V(I)) DIRE 620 DERIV(I+3)=0.5/SQRT(V(J)*V(K)) DIRE 630 DERIV(J)=-V(K)*V(I+3)/(4.*(V(J)*V(K))**1.5) DIRE 640 DERIV(K)=-V(J)*V(I+3)/(4.*(V(J)*V(K))**1.5) DIRE 650 SSA(I+3)=0. DIRE 660 DO 7 J=1,6 DIRE 670 DO 7 K=J,6 DIRE 680 X=1. DIRE 690 IF(J.NE.K)X=2. DIRE 700 7 SSA(I+3)=SSA(I+3)+SM(J,K)*DERIV(J)*DERIV(K)*X DIRE 710 SSA(I+3)=SSA(I+3)/((1.-SA(I+3)**2)*(2.*RAD)**2) DIRE 720 IF(SA(I+3).EQ.0.0) GOTO 25 DIRE 730 SA(I+3)=ATAN(SQRT(1.-SA(I+3)**2)/SA(I+3))/(2.*RAD) DIRE 740 GOTO 30 DIRE 750 25 SA(I+3)=90. DIRE 760 30 CONTINUE DIRE 770 IF(V(I+3).EQ.0.)SA(I+3)=90. DIRE 780 IF(SM(I+3,I+3).EQ.0.)SSA(I+3)=0. DIRE 790 5 CONTINUE DIRE 800 DO 12 I=1,6 DIRE 810 IF(SSA(I).LT.0.)SSA(I)=0. DIRE 820 12 CONTINUE DIRE 830 CALL ESD(SM,V,1.) DIRE 840 WRITE(NWD,8) DIRE 850 DO 9 I=1,3 DIRE 860 I1=10000.*SM(I,I)+0.5 DIRE 870 I2=1000.*SM(I+3,I+3)+0.5 DIRE 880 I3=1000000.*SQRT(SSA(I))+0.5 DIRE 890 I4=1000.*SQRT(SSA(I+3))+0.5 DIRE 900 9 WRITE(NWD,10)V(I),I1,V(I+3),I2,SA(I),I3,SA(I+3),I4 DIRE 910 X1=1./SQRT(VOLUME) DIRE 920 I2=100.*SQRT(SVOL/(4.*VOLUME**3))+0.5 DIRE 930 X3=SQRT(VOLUME) DIRE 940 I4=10000000.*SQRT(SVOL/(4.*VOLUME))+0.5 DIRE 950 WRITE(NWD,11)X1,I2,X3,I4 DIRE 960 8 FORMAT(/14X,6HDIRECT,30X,10HRECIPROCAL/ 7X,4HAXES,13X,6HANGLES,16XDIRE 970 1,4HAXES,14X,6HANGLES) DIRE 980 10 FORMAT(1X,F8.4,1H(,I6 ,1H),F10.3,1H(,I6 ,1H),F13.6,1H(,I6 ,1H),DIRE 990 1F10.3,1H(,I6 ,1H)) DIRE1000 11 FORMAT(/7X,7HVOLUME=,F8.2,1H(,I6 ,1H),14X,7HVOLUME=,F9.7,1H(,I6 DIRE1010 1,1H)) DIRE1020 WRITE(NWD,91) 91 FORMAT(102H0THERMAL B-FACTORS AND STANDARD DEVIATIONS. B(I,J)=BIJ 1/A(I)STAR*A(J)STAR*0.25 = 8*PI*PI*M.S.AMPLITUDE// 2 5H ATOM, 9X,6HB(1,1),13X,6HB(2,2),13X,6HB(3,3),13X,6HB(1,2),13X, 3 6HB(1,3),13X,6HB(2,3)) FAK(1)=1./(SA(1)*SA(1)*0.25) FAK(2)=1./(SA(2)*SA(2)*0.25) FAK(3)=1./(SA(3)*SA(3)*0.25) FAK(4)=1./(SA(1)*SA(2)*0.25) FAK(5)=1./(SA(1)*SA(3)*0.25) FAK(6)=1./(SA(2)*SA(3)*0.25) DO81I=1,N DO82J=1,6 BIJ (J)=XL(I,J+8)*FAK(J) 82 SBIJ(J)=SZ(I,J+8)*FAK(J) 81 WRITE(NWD,92) ATEXT(I),(BIJ(J),SBIJ(J),J=1,6) 92 FORMAT( 1X,A4,3X,F7.4,1H(,F7.4,1H),3X,F7.4,1H(,F7.4,1H),3X,F7.4,1H 1(,F7.4,1H),3X,F7.4,1H(,F7.4,1H),3X,F7.4,1H(,F7.4,1H),3X,F7.4,1H(, 2F7.4,1H)) RETURN DIRE1030 END DIRE1040 SUBROUTINE ESD(SM,V,SUM) ESD( 0 C-----THIS SUBROUTINE CALCULATES THE E.S.D OF THE DIRECT CELL DIMENSIONSESD( 10 C-----FROM THE INVERSE NORMAL MATRIX SM ,THE CELL CONSTANTS V AND THE ESD( 20 C-----REDUCED CHI SQUARED VALUE SUM. ESD( 30 DIMENSION SM(6,6),V(6),S(6),R(6),E(6) ESD( 40 I=0 ESD( 50 DO 1 IA=1,3 ESD( 60 I=I+1 ESD( 70 IB=IA+1 ESD( 80 IF(IB.GT.3)IB=IB-3 ESD( 90 IC=IB+1 ESD( 100 IF(IC.GT.3)IC=IC-3 ESD( 110 ID=IA+3 ESD( 120 IE=IB+3 ESD( 130 IP=IC+3 ESD( 140 FNUM=4.*V(IB)*V(IC)-V(ID)*V(ID) ESD( 150 DEN=4.*V(IA)*V(IB)*V(IC)-V(IA)*V(ID)*V(ID)-V(IB)*V(IE)*V(IE)-V(IC)ESD( 160 1*V(IP)*V(IP)+V(ID)*V(IE)*V(IP) ESD( 170 S(I)=SQRT(FNUM/DEN) ESD( 180 R(IA)=-S(I)**4 ESD( 190 DENL=1./(DEN*DEN) ESD( 200 R(IB)=DENL*(4.*V(IC)*DEN-(4.*V(IA)*V(IC)-V(IE)*V(IE))*FNUM) ESD( 210 R(IC)=DENL*(4.*V(IB)*DEN-(4.*V(IA)*V(IB)-V(IP)*V(IP))*FNUM) ESD( 220 R(ID)=-DENL*(2.*V(ID)*DEN-(2.*V(IA)*V(ID)-V(IE)*V(IP))*FNUM) ESD( 230 R(IE)=DENL*FNUM*(2.*V(IB)*V(IE)-V(ID)*V(IP)) ESD( 240 R(IP)=DENL*FNUM*(2.*V(IC)*V(IP)-V(ID)*V(IE)) ESD( 250 E(I)=ERROR(SM,R,SUM) ESD( 260 E(I)=E(I)/(2.*S(I)) ESD( 270 FNUM=V(IE)*V(IP)-2.*V(IA)*V(ID) ESD( 280 DEN=16.*V(IA)*V(IA)*V(IB)*V(IC)+V(IE)*V(IE)*V(IP)*V(IP)-4.*(V(IA)*ESD( 290 1V(IB)*V(IE)*V(IE)+V(IC)*V(IA)*V(IP)*V(IP)) ESD( 300 DENL=1./(DEN*DEN) ESD( 310 FNUML=FNUM*FNUM ESD( 320 S(I+3)=FNUML/DEN ESD( 330 R(IA)=-4.*FNUM*DENL*(V(ID)*DEN+FNUM*(8.*V(IA)*V(IB)*V(IC)-V(IB)*V(ESD( 340 1IE)*V(IE)-V(IC)*V(IP)*V(IP))) ESD( 350 R(IB)=-4.*FNUML*DENL*(4.*V(IA)*V(IA)*V(IC)-V(IA)*V(IE)*V(IE)) ESD( 360 R(IC)=-4.*FNUML*DENL*(4.*V(IA)*V(IA)*V(IB)-V(IA)*V(IP)*V(IP)) ESD( 370 R(ID)=-4.*FNUM/DEN*V(IA) ESD( 380 R(IE)=2.*FNUM*DENL*(V(IP)*DEN-FNUM*(V(IE)*V(IP)*V(IP)-4.*V(IA)*V(IESD( 390 1B)*V(IE))) ESD( 400 R(IP)=2.*FNUM*DENL*(V(IE)*DEN-FNUM*(V(IE)*V(IE)*V(IP)-4.*V(IC)*V(IESD( 410 1A)*V(IP))) ESD( 420 E(I+3)=ERROR(SM,R,SUM) ESD( 430 IF(S(I+3).EQ.0.0) GOTO 3 ESD( 440 E(I+3)=E(I+3)/(2.*SQRT(S(I+3)*(1.-S(I+3))))*180./3.14159265359 ESD( 450 S(I+3)=ATAN(SQRT((1.-S(I+3))/S(I+3)))*180./3.14159265359 ESD( 460 GOTO 4 ESD( 470 3 S(I+3)=90. ESD( 480 E(I+3)=0.0 ESD( 490 4 CONTINUE ESD( 500 IF(FNUM.LT.0.)S(I+3)=180.-S(I+3) ESD( 510 IF(SM(ID,ID).EQ.0.)E(ID)=0. ESD( 520 1 CONTINUE ESD( 530 IF(ABS(S(1)-S(2)).LT..00008)E(1)=E(2) ESD( 540 DO 2 I=1,6 ESD( 550 SM(I,I)=E(I) ESD( 560 2 V(I)=S(I) ESD( 570 RETURN ESD( 580 END ESD( 590 FUNCTION ERROR(A,B,OMEGA) ERRO 0 C-----THIS FUNCTION CALCULATES THE ESD OF A FUNCTION OF 6 VARIABLES OF ERRO 10 C-----WHICH THE COVARIANCE MATRIX A*OMEGA AND THEIR DERIVATIVES B ARE ERRO 20 C-----GIVEN. ERRO 30 DIMENSION A(6,6),B(6) ERRO 40 SUM=0. ERRO 50 DO 1 I=1,6 ERRO 60 DO 1 J=I,6 ERRO 70 X=2. ERRO 80 IF(I.EQ.J)X=1. ERRO 90 1 SUM=SUM+A(I,J)*B(I)*B(J)*X ERRO 100 IF(SUM.LT.0.)SUM=0. ERRO 110 ERROR=SQRT(SUM*OMEGA) ERRO 120 RETURN ERRO 130 END ERRO 140 SUBROUTINE OUTPUT OUTP 0 C-----THIS SUBROUTINE PRINTS ALL RELEVANT DATA AFTER EACH CYCLE OUTP 10 COMMON CUBIC,FORCED,HALPER,PAC,PUNCH,SKIP,STAR,TYPE,UNIAX, OUTP 20 1VERT,YES OUTP 30 COMMON AH,B1,B2,BH,CH,DELTA,EPS,FACCOS,FACCUB,FACSIN,FLOG,FLOG2, OUTP 40 1FMAGDI,FMAGOB,FNUCDI,FNUCOB,OMEGA,PI,RELAXB,RELAXC,RELAXS,SCALE, OUTP 50 2SKEW,SOM,SUMCAL,SUMDFF,SUMDIF,SUMFOB,SUMOBS,SUMW,SUMWBS,SW,W,XN, OUTP 60 3YCALC,YOBS,ZERO,SLABDA,SLOG,ARG OUTP 70 COMMON YALFA,ICENT,ICYCLE,YOMEG,IORD1,IORD2,IPONS,IPREV,IRL,YSTAP,OUTP 80 1KL,KM,KSTART,LCOORD,LCORD,LIM,LIMO,LLL,LMAT,LOUT,MAXCON,MAXS,MN,N,OUTP 90 2NCYCLE,NKEY,NLC,NN,NQC,NTEXT,NUM,YSS,NO OUTP 100 COMMON ICODE(4889),ICON(4,9),JCON(4,9),LCONI(4,9),LCONJ(4,9),LN(4) 1,LP(71,14),MEQ(69),MTYP(69),NTYP(69),IOBS(10),ICALC(10) COMMON A(71,14),AL(3,3),ATEXT(69),B(6),C(4,98),CF(4),CONA(4,9), 1COSA(3),COSAF(3),COSAR(69,48),COSARI(48),DERIV(149),DERSTO(98,149) 2,DF(4),DUMP(98),DUMPA(500),F(2,21),FAC(4889),FESD(4889),FF(4889,2) 3,FLAG(71,14),FMAG(4889),FNUC(4889),FOBS(4889),H(48,3),HH(4889,3), 4HL(3),HN(3),HW(98),OCC(69),PAK(4889),PH(3),RM(48,9,9),RSA(69,3,3), 5RSB(69,3,3),S(2,21),SA(69),SB(69),SIDE(3),SINAR(69,48),SINARI(48), 6SJJ(4889),SM(48,3,3),SMM(149,149),SMOM(3),SNEX(69),SOME(98), 7SOMEGA(4889),T(48,3),TANN(98),TEMP(69),TEXT(20),THETA(4889),TR(48) 8,V(149),XL(71,14),SZ(71,14) COMMON NRD,NWD,NSD,NID,NFD,NDD,NCD,PHASEN(4889),SP(6,6),RFAC(6) 1,URABS,ABSOR LOGICAL HALPER,PAC,UNIAX,CUBIC,PUNCH,TYPE,STAR,SKIP, OUTP 220 1YES,VERT,FORCED OUTP 230 DIMENSION SY(14) OUTP 240 WRITE(NWD,1) OUTP 250 WRITE(NWD,2)ICYCLE OUTP 260 IF(FORCED)WRITE(NWD,3) OUTP 270 WRITE(NWD,4) OUTP 280 DO 5 I=1,N OUTP 290 DO 65 J=1,4 OUTP 300 FACTOR=RELAXB OUTP 310 IF(J.LT.4)FACTOR=RELAXC OUTP 320 K=LP(I,J) OUTP 330 IF(K)7,8,7 OUTP 340 8 X=0. OUTP 350 Y=0. OUTP 360 Z=0. OUTP 370 GOTO 9 OUTP 380 7 X=A(I,J) OUTP 390 Y=V(K)*X*FACTOR OUTP 400 Z=SQRT(ABS(SMM(K,K))*ABS(X)) ! A.HEWAT 9/6/87 OUTP 410 9 XL(I,J)=XL(I,J)+Y OUTP 420 SY(J)=Y OUTP 430 65 SZ(I,J)=Z OUTP 440 5 WRITE(NWD,6)ATEXT(I),(XL(I,J),SY(J),SZ(I,J),J=1,4) OUTP 450 WRITE(NWD,400) OUTP 455 DO500I=1,N DO650J=9,14 FACTOR=RELAXB K=LP(I,J) IF(K)700,800,700 800 X=0. Y=0. Z=0. GO TO 900 700 X=A(I,J) Y=V(K)*X*FACTOR Z=SQRT(SMM(K,K))*ABS(X) 900 XL(I,J)=XL(I,J)+Y SY(J)=Y 650 SZ(I,J)=Z 500 WRITE(NWD,600) ATEXT(I),(XL(I,J),SY(J),SZ(I,J),J=9,14) OUTP 456 WRITE(NWD,11) OUTP 460 DO 12 I=1,N OUTP 470 DO 13 J=5,8 OUTP 480 K=LP(I,J) OUTP 490 IF(K)14,15,14 OUTP 500 15 X=0. OUTP 510 Y=0. OUTP 520 Z=0. OUTP 530 GOTO 16 OUTP 540 14 X=A(I,J) OUTP 550 Y=V(K)*X*RELAXS OUTP 560 Z=SQRT(SMM(K,K))*ABS(X) OUTP 570 16 XL(I,J)=XL(I,J)+Y OUTP 580 IF((MTYP(I).EQ.0).AND.(J.GT.5))GOTO 17 OUTP 590 SY(J)=Y OUTP 600 SZ(I,J)=Z OUTP 610 GOTO 13 OUTP 620 17 XL(I,J)=0. OUTP 630 SY(J)=0. OUTP 640 SZ(I,J)=0. OUTP 650 13 CONTINUE OUTP 660 IF(MTYP(I).EQ.0)GOTO 18 OUTP 670 Y=0. OUTP 680 Z=0. OUTP 690 DO 21 J=1,3 OUTP 700 K=J+1 OUTP 710 IF(K.GT.3)K=K-3 OUTP 720 L=K+1 OUTP 730 IF(L.GT.3)L=L-3 OUTP 740 21 SMOM(J)=XL(I,J+5)+XL(I,L+5)*COSA(K)+XL(I,K+5)*COSA(L) OUTP 750 DO 22 J=1,3 OUTP 760 DO 22 K=1,3 OUTP 770 LPJ=LP(I,J+5) OUTP 780 LPK=LP(I,K+5) OUTP 790 X1=0. OUTP 800 IF(.NOT.((LP(I,J+5).EQ.0).OR.(LP(I,K+5).EQ.0)))X1=SMM(LPJ,LPK)* OUTP 810 1A(I,J+5)*A(I,K+5) OUTP 820 X2=1. OUTP 830 JK=6-J-K OUTP 840 IF(J.NE.K)X2=COSA(JK) OUTP 850 Y=Y+XL(I,J+5)*XL(I,K+5)*X2 OUTP 860 22 Z=Z+SMOM(J)*SMOM(K)*X1 OUTP 870 IF (Y.EQ.0.0) Y = 1.0E-10 ! J.K.C. Z=SQRT(Z/Y) OUTP 880 Y=SQRT(Y) OUTP 890 YY=Y OUTP 900 ZZ=Z OUTP 910 GOTO 12 OUTP 920 18 YY=0. OUTP 930 ZZ=0. OUTP 940 12 WRITE(NWD,10)ATEXT(I),(XL(I,J),SY(J),SZ(I,J),J=5,8),YY,ZZ OUTP 950 K=LP(N+1,5) OUTP 960 IF(K)25,26,25 OUTP 970 26 X=0. OUTP 980 Y=0. OUTP 990 Z=0. OUTP1000 GOTO 27 OUTP1010 25 X=A(N+1,5) OUTP1020 Y=V(K)*X*RELAXS OUTP1030 Z=SQRT(SMM(K,K))*ABS(X) OUTP1040 27 XL(N+1,5)=XL(N+1,5)+Y OUTP1050 X=1./XL(N+1,5) OUTP1060 Y=X-SCALE OUTP1070 Z=Z/(XL(N+1,5)*XL(N+1,5)) OUTP1080 SCALE=X OUTP1090 WRITE(NWD,28)X,Y,Z OUTP1100 K=LP(N+1,4) OUTP1110 IF(K)29,30,29 OUTP1120 30 X=0. OUTP1130 Y=0. OUTP1140 Z=0. OUTP1150 GOTO 31 OUTP1160 29 X=A(N+1,4) OUTP1170 Y=V(K)*X*RELAXB OUTP1180 Z=SQRT(SMM(K,K))*ABS(X) OUTP1190 31 XL(N+1,4)=XL(N+1,4)+Y OUTP1200 WRITE(NWD,32)XL(N+1,4),Y,Z OUTP1210 WRITE(NWD,33) OUTP1220 DO 34 J=1,3 OUTP1230 K=LP(N+2,J) OUTP1240 IF(K)35,36,35 OUTP1250 36 X=0. OUTP1260 Y=0. OUTP1270 Z=0. OUTP1280 GOTO 37 OUTP1290 35 X=A(N+2,J) OUTP1300 Y=V(K)*X*0.8 OUTP1310 Z=SQRT(SMM(K,K))*ABS(X) OUTP1320 37 XL(N+2,J)=XL(N+2,J)+0.5*Y !A.Hewat 24 May 88 OUTP1330 WRITE(NWD,38)XL(N+2,J),Y,Z OUTP1340 34 CONTINUE OUTP1350 K=LP(N+2,4) OUTP1360 IF(K)39,40,39 OUTP1370 40 X=0. OUTP1380 Y=0. OUTP1390 Z=0. OUTP1400 GOTO 41 OUTP1410 39 X=A(N+2,4) OUTP1420 Y=V(K)*X*0.8 OUTP1430 Z=SQRT(SMM(K,K))*ABS(X) OUTP1440 41 XL(N+2,4)=XL(N+2,4)+0.5*Y !A.Hewat 24 May 88 OUTP1450 WRITE(NWD,42)XL(N+2,4),Y,Z OUTP1460 WRITE(NWD,43) OUTP1470 DO 44 J=1,6 OUTP1480 K=LP(N+2,J+4) OUTP1490 IF(K)45,46,45 OUTP1500 46 X=0. OUTP1510 Y=0. OUTP1520 Z=0. OUTP1530 GOTO 47 OUTP1540 45 X=A(N+2,J+4) OUTP1550 Y=V(K)*X*0.8 OUTP1560 Z=SQRT(SMM(K,K))*ABS(X) OUTP1570 47 XL(N+2,J+4)=XL(N+2,J+4)+Y OUTP1580 44 WRITE(NWD,48)XL(N+2,J+4),Y,Z OUTP1590 K=LP(N+1,1) OUTP1600 IF(K)49,50,49 OUTP1610 50 X=0. OUTP1620 Y=0. OUTP1630 Z=0. OUTP1640 GOTO 51 OUTP1650 49 X=A(N+1,1) OUTP1660 Y=V(K)*X*RELAXB OUTP1670 Z=SQRT(SMM(K,K))*ABS(X) OUTP1680 51 XL(N+1,1)=XL(N+1,1)+Y OUTP1690 WRITE(NWD,52)XL(N+1,1),Y,Z OUTP1700 K=LP(N+1,2) OUTP1710 IF(K)53,54,53 OUTP1720 54 X=0. OUTP1730 Y=0. OUTP1740 Z=0. OUTP1750 GOTO 55 OUTP1760 53 X=A(N+1,2) OUTP1770 Y=V(K)*X*0.8 OUTP1780 Z=SQRT(SMM(K,K))*ABS(X) OUTP1790 55 XL(N+1,2)=XL(N+1,2)+Y OUTP1800 WRITE(NWD,56)XL(N+1,2),Y,Z OUTP1810 RFAC(1)=100.*SUMDFF/SUMFOB OUTP1820 RFAC(2)=100.*SUMDIF/SUMOBS OUTP1830 RFAC(3)=100.*SQRT(SUMW/SUMWBS) OUTP1840 WRITE(NWD,57) (RFAC(I),I=1,3) OUTP1850 WRITE(*,57) (RFAC(I),I=1,3) RFAC(4)=100.*SQRT((FLOAT(NUM)-FLOAT(MAXS)+FLOAT(MAXCON))*250. 1 /SUMWBS) OUTP1860 WRITE(NWD,58) RFAC(4) OUTP1870 RFAC(5)=100.*FNUCDI/FNUCOB OUTP1880 WRITE(NWD,59)RFAC(5) OUTP1890 IF(.NOT.(HALPER.OR.UNIAX.OR.CUBIC))GOTO 60 OUTP1900 IF (FMAGOB.EQ.0.0) FMAGOB = 1.0E-10 ! JKC RFAC(6)=100.*FMAGDI/FMAGOB OUTP1910 WRITE(NWD,61)RFAC(6) OUTP1920 60 I=NUM-MAXS+MAXCON OUTP1930 WRITE(NWD,62)I OUTP1940 X=SUMWBS/250. OUTP1950 Y=SUMW/(250.*FLOAT(I)) OUTP1960 WRITE(NWD,63)SUMDIF,SUMOBS,SUMCAL,X,SUMDFF,SUMFOB,FNUCDI,FNUCOB, OUTP1970 1FMAGDI,FMAGOB,Y,SKEW OUTP1980 1 FORMAT(129H0++++++++++++++++++++++++++++++++++++++++++++++++++++++OUTP1990 1++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++OUTP2000 1++++++++) OUTP2010 2 FORMAT(//,14H CYCLE NUMBER=,I4) OUTP2020 3 FORMAT(21H *FORCED TERMINATION*) OUTP2030 4 FORMAT(/49H NEW PARAMETERS, SHIFTS, AND STANDARD DEVIATIONS=,//107OUTP2040 1H ATOM X DX SX Y DY SY Z OUTP2150 1 DZ SZ B DB SB) OUTP2060 6 FORMAT(1X,A4,3X,4(3F8.5,3X)) OUTP2070 400 FORMAT(129H0ATOM B11 DB11 SB11 B22 DB22 SB22 B33 OUTP2075 1DB33 SB33 B12 DB12 SB12 B13 DB13 SB13 B23 DB2OUTP2076 23 SB23) OUTP2077 600 FORMAT(1X,A4,18F7.4) OUTP2078 10 FORMAT(1X,A4,3X,4(3F8.5,2X),2F8.5) OUTP2080 11 FORMAT( 127H0ATOM N DN SN KX DKX OUTP2090 1SKX KY DKY SKY KZ DKZ SKZ M OUTP2100 1 SM) OUTP2110 28 FORMAT( 22H0OVERALL SCALE FACTOR=,3F8.4) OUTP2120 32 FORMAT(22H OVERALL TEMP. FACTOR=,3F8.4) OUTP2130 33 FORMAT( 22H0HALFWIDTH PARAMETERS=) OUTP2140 38 FORMAT(1X,3F11.2) OUTP2150 42 FORMAT( 11H0ZEROPOINT=,3F7.3) OUTP2160 43 FORMAT( 16H0CELL CONSTANTS=) OUTP2170 48 FORMAT(1X,3F11.7) OUTP2180 52 FORMAT( 33H0PREFERRED ORIENTATION PARAMETER=/3F11.5) OUTP2190 56 FORMAT( 21H0ASYMMETRY PARAMETER=/3F11.5) OUTP2200 57 FORMAT( 11H R-FACTORS=,3F8.2) OUTP2210 58 FORMAT(11H EXPECTED =,8X,F8.2) OUTP2220 59 FORMAT( 20H0R-FACTOR(NUCLEAR) =,F8.2) OUTP2230 61 FORMAT(20H R-FACTOR(MAGNETIC)=,F8.2) OUTP2240 62 FORMAT( 7H0N-P+C=,I8) OUTP2250 63 FORMAT(/,7X,7HSUMYDIF, 7X,7HSUMYOBS, 6X,8HSUMYCALC, 6X,10HSUMWYOBSOUTP2260 1SQ, 5X,7HSUMIDIF, 7X,7HSUMIOBS,/,1X,6E14.4,//,6X,9HSUMNUCDIF, 5X,9OUTP2270 1HSUMNUCOBS, 5X,9HSUMMAGDIF, 5X,9HSUMMAGOBS, 6X,8HRESIDUAL, 6X,8HSKOUTP2280 1EWNESS/1X,6E14.4/) OUTP2290 RETURN OUTP2300 END OUTP2310 SUBROUTINE EXPUT EXPU 0 C-----THIS SUBROUTINE OUTPUTS THE DATA ON PUNCHTAPES EXPU 10 COMMON CUBIC,FORCED,HALPER,PAC,PUNCH,SKIP,STAR,TYPE,UNIAX, EXPU 20 1VERT,YES EXPU 30 COMMON AH,B1,B2,BH,CH,DELTA,EPS,FACCOS,FACCUB,FACSIN,FLOG,FLOG2, EXPU 40 1FMAGDI,FMAGOB,FNUCDI,FNUCOB,OMEGA,PI,RELAXB,RELAXC,RELAXS,SCALE, EXPU 50 2SKEW,SOM,SUMCAL,SUMDFF,SUMDIF,SUMFOB,SUMOBS,SUMW,SUMWBS,SW,W,XN, EXPU 60 3YCALC,YOBS,ZERO,SLABDA,SLOG,ARG EXPU 70 COMMON YALFA,ICENT,ICYCLE,YOMEG,IORD1,IORD2,IPONS,IPREV,IRL,YSTAP,EXPU 80 1KL,KM,KSTART,LCOORD,LCORD,LIM,LIMO,LLL,LMAT,LOUT,MAXCON,MAXS,MN,N,EXPU 90 2NCYCLE,NKEY,NLC,NN,NQC,NTEXT,NUM,YSS,NO EXPU 100 COMMON ICODE(4889),ICON(4,9),JCON(4,9),LCONI(4,9),LCONJ(4,9),LN(4) 1,LP(71,14),MEQ(69),MTYP(69),NTYP(69),IOBS(10),ICALC(10) COMMON A(71,14),AL(3,3),ATEXT(69),B(6),C(4,98),CF(4),CONA(4,9), 1COSA(3),COSAF(3),COSAR(69,48),COSARI(48),DERIV(149),DERSTO(98,149) 2,DF(4),DUMP(98),DUMPA(500),F(2,21),FAC(4889),FESD(4889),FF(4889,2) 3,FLAG(71,14),FMAG(4889),FNUC(4889),FOBS(4889),H(48,3),HH(4889,3), 4HL(3),HN(3),HW(98),OCC(69),PAK(4889),PH(3),RM(48,9,9),RSA(69,3,3), 5RSB(69,3,3),S(2,21),SA(69),SB(69),SIDE(3),SINAR(69,48),SINARI(48), 6SJJ(4889),SM(48,3,3),SMM(149,149),SMOM(3),SNEX(69),SOME(98), 7SOMEGA(4889),T(48,3),TANN(98),TEMP(69),TEXT(20),THETA(4889),TR(48) 8,V(149),XL(71,14),SZ(71,14) COMMON NRD,NWD,NSD,NID,NFD,NDD,NCD,PHASEN(4889),SP(6,6),RFAC(6) 1,URABS,ABSOR LOGICAL HALPER,PAC,UNIAX,CUBIC,PUNCH,TYPE,STAR,SKIP, EXPU 220 1YES,VERT,FORCED EXPU 230 DIMENSION IN(3),SX(10),ISZ(6),ISC(6) EXPU 240 DIMENSION RTENS(6),REELC(6) IF(.NOT.YES)GOTO 1 EXPU 250 WRITE(NID,2)IORD2 EXPU 260 DO 3 I=1,IORD2 EXPU 270 DO30J=1,3 30 IN(J)=HH(I,J)*HN(J) MLT=SJJ(I) X=THETA(I)/100. EXPU 280 3 WRITE(NID,4) ICODE(I),(IN(J),J=1,3),MLT,X,FNUC(I),FMAG(I) EXPU 290 WRITE(NID,5) 1 IF(.NOT.PUNCH)GOTO 6 EXPU 310 C WRITE(NFD,2)IORD2 ! A.Hewat 31Jan90 EXPU 320 DO 50 I=1,3 EXPU 330 50 IN(I)=HN(I) EXPU 340 c WRITE(NFD,7)(IN(I),I=1,3)! A.Hewat 31Jan90 EXPU 350 DO 8 I=1,IORD2 EXPU 360 BB=SJJ(I)*ARG EXPU 370 DO 9 J=1,3 EXPU 380 9 IN(J)=HH(I,J)*HN(J) EXPU 390 C X=FOBS(I)/(FNUC(I)+FMAG(I)) EXPU 40 C...JKC DIV BY ZERO AJJ = FNUC(I)+FMAG(I) IF (AJJ.EQ.0.0) AJJ = 1.0E-10 X=SQRT(ABS(FOBS(I))/AJJ) ! A.Hewat 30Jan90 c X1=FNUC(I)/BB X1=SQRT(4.*FNUC(I)/(BB*IRL*IRL)) ! A.Hewat 30Jan90 X2=X1*X c X3=FMAG(I)/BB X3=SQRT(4.*FMAG(I)/(BB*IRL*IRL)) ! A.Hewat 30Jan90 X4=X3*X C 8 WRITE(NFD,11)(IN(J),J=1,3),X1,X2,X3,X4,PHASEN(I) !A.Hewat 30Jan90 8 WRITE(NFD,11)(IN(J),J=1,3),X1,PHASEN(I),X2,X3,X4 WRITE(NFD,5) 6 IF(LCOORD.NE.1)GOTO 12 WRITE(NDD,13)NCYCLE,RELAXC,RELAXB,RELAXS,LOUT,IPONS,LMAT,LCORD, + LCOORD X=1/XL(N+1,5) EXPU 500 WRITE(NDD,14)X,XL(N+1,4),URABS EXPU 510 FLGD=' ' DO 15 I=1,N EXPU 520 IF(I.EQ.N) FLGD='/' IF(MTYP(I).NE.0)GOTO 15 EXPU 530 DO 115 J=6,8 EXPU 540 115 XL(I,J)=0. EXPU 550 DO116J=1,6 K=LP(N+2,J+4) ISZ(J)=10000.*SZ(I,J)+0.50 ! Hewat 2 Nov 87 116 ISC(J)=10000.*SP(J,J)+0.50 15 WRITE(NDD,16)ATEXT(I),NTYP(I),MTYP(I),MEQ(I),(XL(I,J),J=1,14), EXPU 560 1 (ISZ(J),J=1,5),FLGD WRITE(NDD,18) (XL(N+2,J),J=1,3),XL(N+2,4) 18 FORMAT(3F8.1,F8.4,47X,'/') C.... Always write out real cell parameters!!!! LUN=-1 WRITE(NDD,1718) LUN 1718 FORMAT(I8,71X,'/') DO 1719 J=1,6 RTENS(J)=XL(N+2,J+4) 1719 CONTINUE CALL RCELL(RTENS,REELC) WRITE(NDD,1818) REELC,(ISC(J),J=1,6) 1818 FORMAT(3F8.4,3F8.3,6I4,7X,'/') WRITE(NDD,1918) XL(N+1,1),XL(N+1,2) 1918 FORMAT(2F8.5,63X,'/') WRITE(NDD,2018) MAXS,RFAC 2018 FORMAT(I8,6X,'R-FACTORS=',6F8.4,7X,'/') IF(MAXS.EQ.0)GOTO 19 EXPU 590 FLAGD=' ' DO 20 I=1,N EXPU 600 IF(I.EQ.N) FLAGD='/' IF(MTYP(I).NE.0)GOTO 20 EXPU 610 DO 120 J=6,8 EXPU 620 120 FLAG(I,J)=0. EXPU 630 20 WRITE(NDD,21)(FLAG(I,J),J=1,14),FLAGD EXPU 640 WRITE(NDD,23)FLAG(N+1,5),FLAG(N+1,4),(FLAG(N+2,J),J=1,3),FLAG(N+2,EXPU 650 14),(FLAG(N+2,J),J=5,10),FLAG(N+1,1),FLAG(N+1,2),NLC,NQC EXPU 660 IF(MAXCON.EQ.0)GOTO 19 EXPU 670 DO 24 I=1,MAXCON EXPU 680 FLAGD=' ' K=LN(I) EXPU 690 DO 25 J=1,K EXPU 700 IF(I-NLC)27,27,28 EXPU 710 28 WRITE(NDD,26)CONA(I,J),LCONI(I,J),LCONJ(I,J),ICON(I,J),JCON(I,J) EXPU 720 GOTO 25 EXPU 730 27 WRITE(NDD,26)CONA(I,J),ICON(I,J),JCON(I,J) EXPU 740 25 CONTINUE EXPU 750 CEN=-100. EXPU 760 IF(I.EQ.NLC.OR.I.EQ.MAXCON) FLAGD='/' 24 WRITE(NDD,10) CEN,DF(I),FLAGD EXPU 770 19 WRITE(NDD,5) EXPU 780 12 IF(LCORD.NE.1)GOTO 31 EXPU 790 IL=0 EXPU 800 DO 32 I=1,N EXPU 810 DO 32 J=1,3 EXPU 820 JJ=0 EXPU 830 JJJ=0 EXPU 840 IL=IL+1 EXPU 850 INT=ABS(FLAG(I,J))+0.5 EXPU 860 IF(INT)33,34,33 EXPU 870 34 AA=0 EXPU 880 IB=0 EXPU 890 IC=0 EXPU 900 GOTO 35 EXPU 910 33 AA=A(I,J) EXPU 920 IB=1 EXPU 930 IC=LP(I,J) EXPU 940 35 DO 36 K=1,N EXPU 950 DO 36 L=1,3 EXPU 960 JJ=JJ+1 EXPU 970 INT=ABS(FLAG(K,L))+0.5 EXPU 980 IF(INT)37,38,37 EXPU 990 38 AAL=0. EXPU1000 IBL=0 EXPU1010 ICL=0 EXPU1020 GOTO 39 EXPU1030 37 AAL=A(K,L) EXPU1040 IBL=1 EXPU1050 ICL=LP(K,L) EXPU1060 39 IF(.NOT.((IB.NE.IBL).OR.(IB.EQ.0).OR.(IBL.EQ.0)))GOTO 40 EXPU1070 JJJ=JJJ+1 EXPU1080 SX(JJJ)=0. GOTO 42 40 IF(IL.NE.JJ)GOTO 43 X=AA*AAL*SMM(IC,IC) JJJ=JJJ+1 SX(JJJ)=X GOTO 42 43 X=SMM(ICL,IC) IF(ICL.GE.IC)X=SMM(IC,ICL) X=X*AA*AAL JJJ=JJJ+1 SX(JJJ)=X 42 IF(.NOT.((JJJ.EQ.8).OR.((K.EQ.N).AND.(L.EQ.3)))) GOTO 36 WRITE(NCD,60)(SX(M),M=1,JJJ) JJJ=0 36 CONTINUE 32 CONTINUE 31 WRITE(NWD,5) WRITE(NWD,46) WRITE(NWD,47) 2 FORMAT(I8) 4 FORMAT(5I8,F8.3,2F12.4) 5 FORMAT(/) 7 FORMAT(3I8) 10 FORMAT(F8.1/F8.3) 11 FORMAT(3I5,6F10.4) 13 FORMAT(I8,3F8.4,5I8,7X,1H/) 14 FORMAT(3F8.4,55X,1H/) 16 FORMAT(A4,3I4,8F8.5/6F8.5,5I4,11X,A1) 21 FORMAT(8F8.3/6F8.3,31X,A1) 23 FORMAT(2F8.3,63X,1H//4F8.3,47X,1H//6F8.3,31X,1H//2F8.3,63X, 1 1H//2I8,63X,1H/) 26 FORMAT(F8.3,4I8,39X,A1) 29 FORMAT(2I8) 41 FORMAT(10E9.3) 46 FORMAT(130H0++++++++++++++++++++++++++++++++++++++++++++++++++++++HHHHHHH 1++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++HHHHHHH 2+++++++++) 47 FORMAT(1H1) 60 FORMAT(8E10.4) RETURN END SUBROUTINE CALCUL(NN) C-----THIS SUBROUTINE CALCULATES THE STRUCTURE FACTOR AND ITS DERIVATIVECALC 10 C-----FOR REFLECTION NN ON THE REFLECTION TAPE. CALC 20 COMMON CUBIC,FORCED,HALPER,PAC,PUNCH,SKIP,STAR,TYPE,UNIAX, CALC 30 1VERT,YES CALC 40 COMMON AH,B1,B2,BH,CH,DELTA,EPS,FACCOS,FACCUB,FACSIN,FLOG,FLOG2, CALC 50 1FMAGDI,FMAGOB,FNUCDI,FNUCOB,OMEGA,PI,RELAXB,RELAXC,RELAXS,SCALE, CALC 60 2SKEW,SOM,SUMCAL,SUMDFF,SUMDIF,SUMFOB,SUMOBS,SUMW,SUMWBS,SW,W,XN, CALC 70 3YCALC,YOBS,ZERO,SLABDA,SLOG,ARG CALC 80 COMMON YALFA,ICENT,ICYCLE,YOMEG,IORD1,IORD2,IPONS,IPREV,IRL,YSTAP,CALC 90 1KL,KM,KSTART,LCOORD,LCORD,LIM,LIMO,LLL,LMAT,LOUT,MAXCON,MAXS,MN,N,CALC 100 2NCYCLE,NKEY,NLC,NQ,NQC,NTEXT,NUM,YSS,NO CALC 110 COMMON ICODE(4889),ICON(4,9),JCON(4,9),LCONI(4,9),LCONJ(4,9),LN(4) 1,LP(71,14),MEQ(69),MTYP(69),NTYP(69),IOBS(10),ICALC(10) COMMON A(71,14),AL(3,3),ATEXT(69),B(6),C(4,98),CF(4),CONA(4,9), 1COSA(3),COSAF(3),COSAR(69,48),COSARI(48),DERIV(149),DERSTO(98,149) 2,DF(4),DUMP(98),DUMPA(500),F(2,21),FAC(4889),FESD(4889),FF(4889,2) 3,FLAG(71,14),FMAG(4889),FNUC(4889),FOBS(4889),H(48,3),HH(4889,3), 4HL(3),HN(3),HW(98),OCC(69),PAK(4889),PH(3),RM(48,9,9),RSA(69,3,3), 5RSB(69,3,3),S(2,21),SA(69),SB(69),SIDE(3),SINAR(69,48),SINARI(48), 6SJJ(4889),SM(48,3,3),SMM(149,149),SMOM(3),SNEX(69),SOME(98), 7SOMEGA(4889),T(48,3),TANN(98),TEMP(69),TEXT(20),THETA(4889),TR(48) 8,V(149),XL(71,14),SZ(71,14) COMMON NRD,NWD,NSD,NID,NFD,NDD,NCD,PHASEN(4889),SP(6,6),RFAC(6) 1,URABS,ABSOR LOGICAL HALPER,PAC,UNIAX,CUBIC,PUNCH,TYPE,STAR,SKIP, CALC 230 1YES,VERT,FORCED CALC 240 LOGICAL VERT1 CJH DIMENSION HNN(3),XI(14),AX(3),BX(3),RSAI(3,3),RSBI(3,3),AXI(3),BXICALC 250 1(3),DSA(3),DSB(3),DAX(3),DBX(3) CALC 260 NN10=MOD(NN,LIMO)+1 CALC 270 FNUC(NN)=0. CALC 280 FMAG(NN)=0. CALC 290 PAKNN=0. CALC 300 IF(PAC)PAKNN=PAK(NN) CALC 310 TR(1)=0. CALC 320 DO 1 I=1,3 CALC 330 HNN(I)=HH(NN,I) CALC 340 1 H(1,I)=HNN(I) CALC 350 C-----CALCULATION OF TEMP.FACTOR,POSITION AND HALFWIDTH CALC 360 SS=0. CALC 370 DO 2 I=1,3 CALC 380 DO 2 J=I,3 CALC 390 2 SS=AL(I,J)*HNN(I)*HNN(J)+SS CALC 400 SSNN=0.25*SS CALC 410 TAV=EXP(-2.*(XL(N+1,4)+ABSOR)*SSNN+B1*PAKNN) SINTH=SLABDA*SS CALC 430 COSTH=1.-SINTH CALC 440 TANTH=SQRT(SINTH/COSTH) CALC 450 TANN(NN10)=TANTH CALC 460 THETA(NN)=ATAN(TANTH)*PI/SLABDA CALC 470 BB=SQRT(AH*TANTH*TANTH+BH*TANTH+CH)*SLOG CALC 480 HW(NN10)=BB CALC 490 BB=BB*BB CALC 500 C-----VERT=.TRUE. IF ASYMMETRY CORRECTION IS TO BE CALCULATED CALC 510 VERT=.FALSE. CALC 520 IF(THETA(NN).LE.FLOAT(LIM))VERT=.TRUE. CALC 530 VERT1=LIM.LT.0 CJH IF(IRL-2)3,4,4 CALC 540 4 DO 5 IR=2,IRL CALC 550 X=0. CALC 560 DO 6 I=1,3 CALC 570 Y=0. CALC 580 X=T(IR,I)*HNN(I)+X CALC 590 DO 7 J=1,3 CALC 600 7 Y=SM(IR,J,I)*HNN(J)+Y CALC 610 6 H(IR,I)=Y CALC 620 TR(IR)=X CALC 630 5 CONTINUE CALC 640 C-----CALCULATION OF COS(H.X),SIN(H.X) AND TEMP. FACTOR FOR EACH ATOM CALC 650 3 DO 8 I=1,N CALC 660 DO 9 J=1,14 CALC 670 9 XI(J)=XL(I,J) CALC 680 DO 10 IR=1,IRL CALC 690 ARG=TR(IR) CALC 700 DO 11 J=1,3 CALC 710 11 ARG=H(IR,J)*XI(J)+ARG CALC 720 ARG=6.28318530718*ARG CALC 730 ARG2= H(IR,1)*H(IR, 1)*XI( 9)+ H(IR,2)*H(IR,2)*XI(10) CALC 735 1 + H(IR,3)*H(IR, 3)*XI(11)+2.*H(IR,1)*H(IR,2)*XI(12) CALC 736 2 +2.*H(IR,1)*H(IR, 3)*XI(13)+2.*H(IR,2)*H(IR,3)*XI(14) CALC 737 EXPARG=EXP(-ARG2) CALC 735 COSAR(I,IR)=COS(ARG)*EXPARG CALC 740 10 SINAR(I,IR)=SIN(ARG)*EXPARG CALC 750 TEMP(I)=EXP(-XL(I,4)*SSNN) CALC 760 8 CONTINUE CALC 770 C-----ZEROIZE THE DERIVATIVES OF THIS REFLECTION W.R.T. TO PARAMETERS CALC 780 DO 12 I=1,MAXS CALC 790 12 DERIV(I)=0. CALC 800 IF(ICODE(NN).EQ.2)GOTO 13 CALC 810 C-----CALCULATE NUCLEAR CONTRIBUTION TO S.F. AND DERIVATIVES CALC 820 AV=0. CALC 830 BV=0. CALC 840 DO 14 I=1,N CALC 850 IF(NTYP(I).EQ.0)GOTO 14 CALC 860 SAI=0. CALC 870 SBI=0. CALC 880 DO 15 IR=1,IRL CALC 890 SAI=SAI+COSAR(I,IR) CALC 900 IF(ICENT.EQ.1)SBI=SINAR(I,IR)+SBI CALC 910 15 CONTINUE CALC 920 SAI=FLOAT(ICENT)*SAI CALC 930 SA(I)=SAI CALC 940 SB(I)=SBI CALC 950 NI=NTYP(I) CALC 960 SNEXI=B(NI)*XL(I,5)*TEMP(I) CALC 970 C-----CALCULATE A AND B OF F CALC 980 AV=SNEXI*SAI+AV CALC 990 BV=SNEXI*SBI+BV CALC1000 SNEX(I)=2.*SNEXI*TAV CALC1010 14 CONTINUE CALC1020 FNN=(AV*AV+BV*BV)*TAV CALC1030 FNUC(NN)=FNN CALC1040 IF(ICYCLE.LT.NCYCLE) GO TO 140 PHASEN(NN)=3.14159265 AH258/85 IF(AV.NE.0.) PHASEN(NN)=ATAN(BV/AV) AH258/85 IF(AV.LT.0.) PHASEN(NN)=PHASEN(NN)+3.14159265 AH258/85 IF(PHASEN(NN).LT.0.) PHASEN(NN)=PHASEN(NN)+6.28318530 AH258/85 140 IF(MAXS.EQ.0)GOTO 16 CALC1050 C-----CALCULATE DERIVATIVES CALC1060 N1=N+1 CALC1070 DO 17 I=1,N1 CALC1080 IF(I-N1)18,19,19 CALC1090 19 SKIP=.FALSE. CALC1100 GOTO 20 CALC1110 18 IF(NTYP(I).EQ.0)GOTO 17 CALC1120 SKIP=.TRUE. CALC1130 SNEXI=SNEX(I) CALC1140 SAI=SA(I) CALC1150 SBI=SB(I) CALC1160 DO 21 IR=1,IRL CALC1170 SINARI(IR)=SINAR(I,IR) CALC1180 21 COSARI(IR)=COSAR(I,IR) CALC1190 20 DO 22 J=1,14 CALC1200 K=LP(I,J) CALC1210 IF(K.EQ.0)GOTO 22 CALC1220 IF(J.GT.5) GO TO 221 IF(SKIP)GOTO 23 CALC1230 IF(J-1)24,25,24 CALC1240 25 DER=PAKNN*FNN CALC1250 GOTO 26 CALC1260 24 IF(J.EQ.2)GOTO 22 CALC1270 23 IF(J.GT.3)GOTO 29 CALC1280 SUMA=0. CALC1290 SUMB=0. CALC1300 DO 27 IR=1,IRL CALC1310 SUMA=H(IR,J)*SINARI(IR)+SUMA CALC1320 IF(ICENT.EQ.1)SUMB=H(IR,J)*COSARI(IR)+SUMB CALC1330 27 CONTINUE CALC1340 DER=-(AV*SUMA-BV*SUMB)*SNEXI*FLOAT(ICENT)*6.28318530718 CALC1350 GOTO 26 CALC1360 29 IF(J.GT.4)GOTO 28 CALC1370 IF(SKIP)GOTO 30 CALC1380 DER=-2.*SSNN*FNN CALC1390 GOTO 26 CALC1400 30 DER=-(SAI*AV+SBI*BV)*SSNN*SNEXI CALC1410 GOTO 26 CALC1420 28 IF(SKIP)GOTO 31 CALC1430 DER=SCALE*FNN CALC1440 GOTO 26 CALC1450 31 DER=(SAI*AV+SBI*BV)*SNEXI/XL(I,5) CALC1460 GO TO 26 221 IF(J.LT.9) GO TO 22 IF(.NOT.SKIP) GO TO 22 IF(J.GT.11) GO TO 222 K1=J-8 K2=K1 GO TO 224 222 K1=1 K2=2 IF(J.GE.13) K2=3 IF(J.EQ.14) K1=2 224 SUMA=0. SUMB=0. DO225IR=1,IRL SUMA=SUMA+H(IR,K1)*H(IR,K2)*COSARI(IR) IF(ICENT.EQ.1) SUMB=SUMB+H(IR,K1)*H(IR,K2)*SINARI(IR) 225 CONTINUE DER=-(AV*SUMA+BV*SUMB)*SNEXI*FLOAT(ICENT) IF(K1.NE.K2) DER=2.*DER 26 DERIV(K)=A(I,J)*DER+DERIV(K) CALC1470 22 CONTINUE CALC1480 17 CONTINUE CALC1490 16 IF(ICODE(NN).EQ.1)GOTO 32 CALC1500 C-----CALCULATE MAGNETIC CONTRIBUTION TO S.F. AND DERIVATIVES CALC1510 C-----FOR UNIAXIAL SPIN CONFIGURATION CALC1520 13 IF(.NOT.UNIAX)GOTO 33 CALC1530 X=FAC(NN) CALC1540 FACCOS=0.5*(1.+X) CALC1550 FACSIN=1.-X CALC1560 33 DO 34 I=1,3 CALC1570 AX(I)=0. CALC1580 34 BX(I)=0. CALC1590 DO 35 I=1,N CALC1600 IF(MTYP(I).EQ.0)GOTO 35 CALC1610 MEQI=MEQ(I) CALC1620 MTYPI=MTYP(I) CALC1630 SNEXI=FF(NN,MTYPI)*TEMP(I)*0.2695*XL(I,5) CALC1640 SNEX(I)=SNEXI CALC1650 DO 36 IR=1,IRL CALC1660 SINARI(IR)=SINAR(I,IR) CALC1670 36 COSARI(IR)=COSAR(I,IR) CALC1680 DO 37 J=1,3 CALC1690 DO 37 K=1,3 CALC1700 JK=(J-1)*3+K CALC1710 Y=0. CALC1720 Z=0. CALC1730 DO 38 IR=1,IRL CALC1740 X=RM(IR,MEQI,JK) CALC1750 Y=COSARI(IR)*X+Y CALC1760 38 IF(ICENT.EQ.1)Z=SINARI(IR)*X+Z CALC1770 RSAI(J,K)=Y CALC1780 37 RSBI(J,K)=Z CALC1790 DO 39 J=1,3 CALC1800 DO 39 K=1,3 CALC1810 RSAI(J,K)=RSAI(J,K)*FLOAT(ICENT) CALC1820 RSA(I,J,K)=RSAI(J,K) CALC1830 39 RSB(I,J,K)=RSBI(J,K) CALC1840 DO 40 J=1,3 CALC1850 AXI(J)=0. CALC1860 BXI(J)=0. CALC1870 DO 40 K=1,3 CALC1880 AXI(J)=XL(I,K+5)*RSAI(J,K)+AXI(J) CALC1890 IF(ICENT.EQ.1)BXI(J)=XL(I,K+5)*RSBI(J,K)+BXI(J) CALC1900 40 CONTINUE CALC1910 DO 41 J=1,3 CALC1920 AX(J)=AXI(J)*SNEXI+AX(J) CALC1930 IF(ICENT.EQ.1)BX(J)=BXI(J)*SNEXI+BX(J) CALC1940 41 CONTINUE CALC1950 SNEX(I)=SNEXI*TAV CALC1960 35 CONTINUE CALC1970 C-----ACCORDING TO HALPERN AND JOHNSON CALC1980 IF(.NOT.HALPER)GOTO 42 CALC1990 X=0. CALC2000 Y=0. CALC2010 DO 43 I=1,3 CALC2020 DO 43 J=I,3 CALC2030 K=6-I-J CALC2040 X1=1. CALC2050 IF(I.NE.J)X1=2.*COSA(K) CALC2060 X=X1*AX(I)*AX(J)+X CALC2070 X1=2. CALC2080 IF(I.EQ.J)X1=1. CALC2090 Y=X1*HNN(I)*HNN(J)*AX(I)*AX(J)/(SIDE(I)*SIDE(J))+Y CALC2100 IF(ICENT.NE.1)GOTO 43 CALC2110 X1=1. CALC2120 IF(I.NE.J)X1=2.*COSA(K) CALC2130 X=X1*BX(I)*BX(J)+X CALC2140 X1=2. CALC2150 IF(I.EQ.J)X1=1. CALC2160 Y=X1*HNN(I)*HNN(J)*BX(I)*BX(J)/(SIDE(I)*SIDE(J))+Y CALC2170 43 CONTINUE CALC2180 FNN=X-0.25*Y/SSNN CALC2190 C-----FOR UNIAXIAL SPIN CONFIGURATION CALC2200 42 IF(.NOT.UNIAX)GOTO 44 CALC2210 X=0. CALC2220 Y=2.*COSA(3) CALC2230 DO 45 I=1,2 CALC2240 DO 45 J=I,2 CALC2250 X1=Y CALC2260 IF(I.EQ.J)X1=1. CALC2270 X=X1*AX(I)*AX(J)+X CALC2280 IF(ICENT.EQ.1)X=X1*BX(I)*BX(J)+X CALC2290 45 CONTINUE CALC2300 X1=0. CALC2310 IF(ICENT.EQ.1)X1=BX(3)*BX(3) CALC2320 FNN=(X1+AX(3)*AX(3))*FACSIN+FACCOS*X CALC2330 C-----FOR CUBIC SPIN CONFIGURATION CALC2340 44 IF(.NOT.CUBIC)GOTO 46 CALC2350 X=0. CALC2360 DO 47 I=1,3 CALC2370 X=AX(I)*AX(I)+X CALC2380 IF(ICENT.EQ.1)X=BX(I)*BX(I)+X CALC2390 47 CONTINUE CALC2400 FNN=FACCUB*X CALC2410 C-----CALCULATE MAGNETIC S.F. CALC2420 46 FNN=FNN*TAV CALC2430 C-----TOTAL S.F. CALC2440 FMAG(NN)=FNN CALC2450 IF(MAXS.EQ.0)GOTO 48 CALC2460 C-----CALCULATE DERIVATIVES CALC2470 DO 49 I=1,3 CALC2480 DSA(I)=0. CALC2490 49 DSB(I)=0. CALC2500 IF(.NOT.HALPER)GOTO 50 CALC2510 DO 51 I=1,3 CALC2520 J=I+1 CALC2530 IF(J.GT.3)J=J-3 CALC2540 K=J+1 CALC2550 IF(K.GT.3)K=K-3 CALC2560 X=HNN(I)/SIDE(I) CALC2570 Y=HNN(J)/SIDE(J) CALC2580 Z=HNN(K)/SIDE(K) CALC2590 DSA(I)=2.*(AX(I)+AX(J)*COSA(K)+AX(K)*COSA(J)-0.25/SSNN*X*(X*AX(I)+CALC2600 1Y*AX(J)+Z*AX(K))) CALC2610 IF(ICENT.EQ.1)DSB(I)=2.*(BX(I)+BX(J)*COSA(K)+BX(K)*COSA(J)-0.25/SSCALC2620 1NN*X*(X*BX(I)+Y*BX(J)+Z*BX(K))) CALC2630 51 CONTINUE CALC2640 50 IF(.NOT.UNIAX)GOTO 52 CALC2650 DO 53 I=1,2 CALC2660 J=I+1 CALC2670 IF(J.GT.2)J=J-2 CALC2680 DSA(I)=(AX(J)*COSA(3)+AX(I))*2.*FACCOS CALC2690 IF(ICENT.EQ.1)DSB(I)=(BX(J)*COSA(3)+BX(I))*2.*FACCOS CALC2700 53 CONTINUE CALC2710 DSA(3)=AX(3)*2.*FACSIN CALC2720 IF(ICENT.EQ.1)DSB(3)=BX(3)*2.*FACSIN CALC2730 52 IF(.NOT.CUBIC)GOTO 54 CALC2740 DO 55 I=1,3 CALC2750 DSA(I)=AX(I)*2.*FACCUB CALC2760 IF(ICENT.EQ.1)DSB(I)=BX(I)*2.*FACCUB CALC2770 55 CONTINUE CALC2780 54 N1=N+1 CALC2790 DO 56 I=1,N1 CALC2800 IF(I-N1)57,58,58 CALC2810 58 SKIP=.FALSE. CALC2820 GOTO 59 CALC2830 57 IF(MTYP(I).EQ.0)GOTO 56 CALC2840 SKIP=.TRUE. CALC2850 SNEXI=SNEX(I) CALC2860 MEQI=MEQ(I) CALC2870 DO 60 IR=1,IRL CALC2880 SINARI(IR)=SINAR(I,IR) CALC2890 60 COSARI(IR)=COSAR(I,IR) CALC2900 59 DO 61 J=1,8 CALC2910 IF(LP(I,J).EQ.0)GOTO 61 CALC2920 IF(SKIP)GOTO 62 CALC2930 IF(J-1)63,64,63 CALC2940 64 DER=PAKNN*FNN CALC2950 GOTO 65 CALC2960 63 IF(J.EQ.2)GOTO 61 CALC2970 62 DO 66 K=1,3 CALC2980 DAX(K)=0. CALC2990 66 DBX(K)=0. CALC3000 IF(J.GT.3)GOTO 67 CALC3010 DO 68 L=1,3 CALC3020 DO 69 M=1,3 CALC3030 X=0. CALC3040 Y=0. CALC3050 DO 70 IR=1,IRL CALC3060 LM=(L-1)*3+M CALC3070 Z=RM(IR,MEQI,LM) CALC3080 X=H(IR,J)*SINARI(IR)*Z+X CALC3090 IF(ICENT.EQ.1)Y=H(IR,J)*COSARI(IR)*Z+Y CALC3100 70 CONTINUE CALC3110 DAX(L)=XL(I,M+5)*X+DAX(L) CALC3120 IF(ICENT.EQ.1)DBX(L)=XL(I,M+5)*Y+DBX(L) CALC3130 69 CONTINUE CALC3140 DAX(L)=-DAX(L)*SNEXI*FLOAT(ICENT)*6.28318530718 CALC3150 IF(ICENT.EQ.1)DBX(L)=DBX(L)*SNEXI*6.28318530718 CALC3160 68 CONTINUE CALC3170 GOTO 71 CALC3180 67 IF(J.GT.5)GOTO 72 CALC3190 IF(.NOT.SKIP)GOTO 73 CALC3200 DO 74 K=1,3 CALC3210 DO 75 L=1,3 CALC3220 DAX(K)=RSA(I,K,L)*XL(I,L+5)+DAX(K) CALC3230 IF(ICENT.EQ.1)DBX(K)=RSB(I,K,L)*XL(I,L+5)+DBX(K) CALC3240 75 CONTINUE CALC3250 X1=-SSNN CALC3260 IF(J.EQ.4)X1=-2.*SSNN CALC3330 IF(J.NE.4)X1=1./XL(I,5) CALC3270 DAX(K)=DAX(K)*SNEXI*X1 CALC3280 IF(ICENT.EQ.1)DBX(K)=DBX(K)*SNEXI*X1 CALC3290 74 CONTINUE CALC3300 GOTO 71 CALC3310 73 X1=SCALE CALC3320 DER=FNN*X1 CALC3340 GOTO 65 CALC3350 72 DO 76 K=1,3 CALC3360 DAX(K)=RSA(I,K,J-5)*SNEXI CALC3370 IF(ICENT.EQ.1)DBX(K)=RSB(I,K,J-5)*SNEXI CALC3380 76 CONTINUE CALC3390 71 DER=0. CALC3400 DO 77 K=1,3 CALC3410 X1=0. CALC3420 IF(ICENT.EQ.1)X1=DSB(K)*DBX(K) CALC3430 77 DER=X1+DSA(K)*DAX(K)+DER CALC3440 65 LPIJ=LP(I,J) CALC3450 DERIV(LPIJ)=A(I,J)*DER+DERIV(LPIJ) CALC3460 61 CONTINUE CALC3470 56 CONTINUE CALC3480 32 IF(MAXS.EQ.0)GOTO 48 CALC3490 C-----CALCULATE DERIVATIVES CALC3500 FNN=FNUC(NN)+FMAG(NN) CALC3510 SINTH=FNN*PI/(SQRT(SINTH*COSTH)*BB) CALC3520 SS=FNN/(FLOG2*BB) CALC3530 X=TANTH*TANTH CALC3540 DO 78 J=1,3 CALC3550 K=LP(N+2,J) CALC3560 IF(K.EQ.0)GOTO 78 CALC3570 DERIV(K)=X*A(N+2,J)*SS+DERIV(K) CALC3580 78 X=X/TANTH CALC3590 K=LP(N+2,4) CALC3600 IF(K.NE.0)DERIV(K)=A(N+2,4)*2.*FNN/BB+DERIV(K) CALC3610 DO 79 J=1,6 CALC3620 K=LP(N+2,J+4) CALC3630 IF(K.EQ.0)GOTO 79 CALC3640 IF(J.LT.4)X=HNN(J)*HNN(J) CALC3650 IF(J.EQ.4)X=HNN(2)*HNN(3) CALC3660 IF(J.EQ.5)X=HNN(1)*HNN(3) CALC3670 IF(J.EQ.6)X=HNN(1)*HNN(2) CALC3680 DERIV(K)=X*SINTH*A(N+2,J+4)+DERIV(K) CALC3690 79 CONTINUE CALC3700 K=LP(N+1,2) CALC3710 IF((K.NE.0).AND.VERT)DERIV(K)=-A(N+1,2)*FNN/TANTH+DERIV(K) CALC3720 IF((K.NE.0).AND.VERT1) CJH 1 DERIV(K)=-A(N+1,2)*FNN*(1.-TANTH**2)/(2.*TANTH)+DERIV(K) CJH C-----STORE DERIVATIVES FOR LIMO REFLECTIONS AT A TIME CALC3730 DO 80 I=1,MAXS CALC3740 80 DERSTO(NN10,I)=DERIV(I) CALC3750 48 RETURN CALC3760 END CALC3770 SUBROUTINE TIMER(SECONDS) C ******* C Get the Vax cpu time and local time. Initialise with SECONDS=0 C Include Job Process Information DEFinition codes (JPIDEF). C Supplied by DIGITAL. INCLUDE '($JPIDEF)' INTEGER*4 CLOCK_TICKS,CLOCK_TICKS2,SECONDS LOGICAL*4 STATUS COMMON/CPUTIM/START_TIME,CLOCK_TICKS IF(SECONDS.NE.0) GO TO 2 C Get local time. START_TIME = SECNDS(0.0) C Get cpu time. STATUS = LIB$GETJPI(JPI$_CPUTIM,,,CLOCK_TICKS,,) C Return to calling subroutine. RETURN C Get elapsed time. 2 ETIME = SECNDS(START_TIME) C Get cpu time. STATUS = LIB$GETJPI(JPI$_CPUTIM,,,CLOCK_TICKS2,,) C Total cpu time used. CTIME = FLOAT(CLOCK_TICKS2 - CLOCK_TICKS) / 100.0 TYPE 9, CTIME,ETIME 9 FORMAT(/1X,'CPU time :',F9.2,5X,'Elapsed time :',F9.2/) C Return to the monitor. CALL EXIT C ******* END