C**************************** MAIN PROGRAM *****************************
C MAXIMUM 3000 POINTS IN THE SCAN
C MAXIMUM 41 PARAMETERS VARIABLE  THEREFORE AM(41*42/2)
C MAXIMUM 400 PARAMETERS
C MAXIMUM 40 ATOMS
C MAXIMUM 1000 REFLECTIONS
C MAXIMUM 31 POINTS PER REFLECTION, HALFWIDTH 15
C BLOCKSIZE 2 OF 50
      COMMON IOI,IOO,AM(861),V(41),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB,LATTIS,MINRH,MAXRH,NBRAGG,INDEX(5,1500),STEP,
     2 WISE(3100,2),MINPT,MAXPT,PTMIN,PTMAX,RMAX,ISCREW(3),IGLIDE(3),
     3 NREFS,IC,NC,NV,NVM,NP,NM,ITF,ISTOP,KEY,NPC,NA,LVAR(41),TWOPIS,
     4 BETA(6,40),P(400),KI(400),XYZ(3,40),SF(40),AI(40),TO,PWISE(9),
     5 BH,DADAI(40),DADX(3,40),DADB(6,40),H(3),NF,ICENT,NS,NREFL,RHO,
     6 NOUT,NPM,TWOPI,YC,SINTHL,DC(400),FX(32,10),TS(3,24),IS(2,3,24),
     7 NO,NATOM,NNA,NNB,NNC,NND,NNE,NNF,NNG,NNH,LUM,NVL,YSUM,YMASK,
     8 YAV,NFRAC,NMAX,NEXT,IOF,IW,INHKL
      DIMENSION DWISE(100,41),DIFFAC(6),INDIC(400)
   41 FORMAT(1H047X'NUMERATOR'9X'DENOMINATOR'5X'R')
   42 FORMAT(25H0R FACTOR INCLUDING ZEROS,17X,F15.3,4X,F15.3,4X,F5.3)
   43 FORMAT(3I4,3I6,F12.2,F8.2)
   44 FORMAT(10F12.1)
   45 FORMAT(3I8,2F8.2,I8,8A4)
   46 FORMAT('0AGREEMENT FACTORS BASED ON PARAMETERS BEFORE CYCLE'I3/'0S
     1UM(W*(O-C)**2) IS'E11.3/'0SQRT(SUM(W*(0-C)**2)/(NO-NV)) IS'F10.4)
   47 FORMAT('0ESTIMATED FACTORS BASED ON PARAMETERS AFTER  CYCLE'I3/'0S
     1UM(W*(O-C)**2) IS'E11.3/'0SQRT(SUM(W*(0-C)**2)/(NO-NV)) IS'F10.4)
   48 FORMAT(62H MATRIX HAS A ZERO DIAGONAL ELEMENT CORRESPONDING TO PAR
     1AMETERI3,16H OF THOSE VARIED)
   49 FORMAT(//20(' + + +')/'PARAMETERS AFTER LEAST SQUARES CYCLE'I3/'0
     1PARAMETER'10X'OLD'8X'CHANGE'9X'NEW'9X'ERROR'7X'RATIO')
   50 FORMAT(8X,I6,6F13.7)
   51 FORMAT(19H0CORRELATION MATRIX)
   52 FORMAT(1H0I3,10F9.4/(1H 3X,10F9.4))
C     FORMATS 53, 54 AND 55 ARE WITHIN THE LINAC REPLACEMENT REGION
   56 FORMAT(' LIMITING SPHERE REACHED',I6,' REFLECTIONS USED')
   57 FORMAT(' MEAN F.W.H.H. ='I4/' EQUIVALENT NUMBER OF REFLECTIONS ='I
     14/' NEW SQSIG ='F10.4)
      IOI = 5
      IOO = 6
      IOF = 9
C*********ALTERATION POINT 3 ---- INPUT=2 LABELLER *********************
      I=0
      CALL CONSTR(2,I)
      NVM=41
      NPM=400
      NATOM=40
      NREFS=1500
      NWIDE=15
      KBL=50
      KBLKBL=2*KBL
      NMAX=3100-KBLKBL
      KTEST=KBLKBL-NWIDE
      PI=3.1415926
      TWOPIS=2.0*PI**2
      TWOPI=2.0*PI
      DEGRAD=PI/180.
      GAUSS=2.7726
      CONMEW=1.253314137
      YMASK=-4095.
      INDIC(1)=1
      DO 1001 I = 2,NPM
      K=2
      IF(I.GT.5)  K=3
      IF(I.GT.9)  K=2
      IF(I.GT.15) K=6
 1001 INDIC(I)=K
      INDIC(9)=1
      READ(IOI,45)  MINPT,MAXPT,NBRAGG,STEP,WAVETC,NEXT,(P(I),I=1,8)
      WRITE(IOO,45) MINPT,MAXPT,NBRAGG,STEP,WAVETC,NEXT,(P(I),I=1,8)
C.....LINAC ALTERATION STARTING POINT
      PTMIN=FLOAT(MINPT)*STEP
      PTMAX=FLOAT(MAXPT)*STEP
C     THE LIMITS ARE SET 10 PERCENT TOO WIDE TO ALLOW FOR REFINEMENT
      MINRH= 9000.*(SIN(PTMIN*DEGRAD*0.5)/WAVETC)**2
      RMAX=         SIN(PTMAX*DEGRAD*0.5)/WAVETC
      MAXRH=10600.*RMAX**2
   54 FORMAT(//19X'10000'/'   H   K   L  MULT  *RHO POINT 2*THETA')
   55 FORMAT(' MIN. AND MAX. 2-THETA AND POINT',2F8.2,2I8,/' WAVELENGTH'
     1,F12.4)
      CALL PRELIM
      WRITE(IOO,55) PTMIN,PTMAX,MINPT,MAXPT,WAVETC
      IF(INHKL.EQ.1) GO TO 1004
      WRITE(IOO,54)
      DO 1003 I=1,NBRAGG
      RHO=WAVETC*SQRT(0.0001*INDEX(5,I))
      IF(RHO.LT.0.995) GO TO 1002
      NBRAGG=I-1
      WRITE(IOO,56) NBRAGG
      GO TO 1004
 1002 TTHETA=ARSIN(RHO)*2.0/DEGRAD
      IPOINT=TTHETA/STEP
 1003 WRITE(IOO,43) (INDEX(J,I),J=1,5),IPOINT,TTHETA,RHO
 1004 IF(NREFS.EQ.-1) CALL FINISH(6)
      IF(NC.EQ.0) NC=1
      DO 1036 IC = 1,NC
      DO 1005 I=1,NM
 1005 AM(I)=0.0
      DO 1007 I=1,NVL
      DO 1006 J = 1,KBLKBL
 1006 DWISE(J,I)=0.0
 1007 V(I)=0.0
      DO 1008 I = 1,NMAX
 1008 WISE(I,2)=0.0
      MEANWD=0
      NBLOCK=0
      BLOCK=0.0
      RZNUM=0.0
      SIG=0.0
      IQUERY=0
C     START LOOP THROUGH OBSERVATIONS
      DO 1027 NREFL = 1,NBRAGG
C     ENTER SUBROUTINE CALC TO COMPUTE Y(CALC) AND DERIVATIVES
      CALL CALC
      SINTH=SINTHL*WAVETC
      POINT=2.0*ARSIN(SINTH)/DEGRAD
      COSTH=SQRT(1.0-SINTH**2)
      TANTH=SINTH/COSTH
      UVW=(PWISE(6)*TANTH+PWISE(7))*TANTH+PWISE(8)
      FACTOR=WAVETC/(SINTHL*COSTH*DEGRAD)
      AREA=SQRT(UVW)
C     LORENTZ CORRECTION, POLARISATION MUST BE ADDED FOR X-RAYS
      CORLOR=1./(COSTH*SINTH**2)
      YC=YC*CORLOR
      IF(IC.EQ.NOUT) WRITE(IOO,43) (INDEX(L,NREFL),L=1,5),NREFL,YC,POINT
 1009 PHICAL=POINT-PWISE(5)-BLOCK
      MIDREF=PHICAL/STEP
      INDEX(5,NREFL)=NBLOCK+MIDREF
      IF(MIDREF.GT.KTEST) GO TO 1014
      IQUERY=NREFL
      NWIDTH=2.0*AREA/STEP
      MEANWD=MEANWD+NWIDTH
      IF(NWIDTH.GT.NWIDE) NWIDTH=NWIDE
      NSTART=MIDREF-NWIDTH
      NFINIS=MIDREF+NWIDTH
      IF(NSTART.LT.1) CALL FINISH(2)
      DC(5)=-1.0
      DC(6)= TANTH**2
      DC(7)= TANTH
      DC(8)= 1.0
      DC(9)= 1.0
      DO 1010 L = 1,6
 1010 DC(L+9)=(DCELL(L,1)*H(1)**2+DCELL(L,4)*H(2)*H(3)
     1        +DCELL(L,2)*H(2)**2+DCELL(L,5)*H(3)*H(1)
     2        +DCELL(L,3)*H(3)**2+DCELL(L,6)*H(1)*H(2))*FACTOR
      IF(LATTIS.LT.8) GO TO 1011
      DC(10)=DC(10)+DC(11)
      IF(LATTIS.LT.12) GO TO 1011
      IF(LATTIS.GT.19) GO TO 1011
      DC(10)=DC(10)+DC(12)
      DC(13)=DC(13)+DC(14)+DC(15)
 1011 DO 1013 NSTEP = NSTART,NFINIS
      IBLOCK=NSTEP+NBLOCK
      DELTA=NSTEP*STEP-PHICAL
      DELDEL=DELTA*DELTA
      SKEW=DELDEL/TANTH
      IF(NSTEP.GT.MIDREF) SKEW=-SKEW
      EXPON=GAUSS*DELDEL/UVW
      EXPART=EXP(-EXPON)*INDEX(4,NREFL)*(1.0+PWISE(9)*SKEW)/AREA
      YCPART=YC*EXPART
      DIFFAC(2)=2.0*GAUSS*DELTA*YCPART/UVW
      DIFFAC(3)=YCPART*(EXPON-0.5)/UVW
      DIFFAC(1)=YCPART*SKEW
      DIFFAC(6)=EXPART*CORLOR
      IF(KEY.GT.NV) GO TO 1013
      DO 1012 L = KEY,NV
      M=LVAR(L)
      K=INDIC(M)
 1012 DWISE(NSTEP,L)=DWISE(NSTEP,L)+DIFFAC(K)*DC(M)
 1013 WISE(IBLOCK,2)=WISE(IBLOCK,2)+YCPART
      IF(IQUERY.LT.NBRAGG) GO TO 1027
C     ADD DIFFERENTIALS FOR BACKGROUND AND SCALE IF REQUIRED AND
C     ADD BLOCK TO LEAST SQUARES MATRIX, SUM DIFFERENCES
 1014 DO 1025 KBLOCK = 1,KBL
      IBLOCK=KBLOCK+NBLOCK
      BACKGD=PWISE(2)
      L=0
      IF(KI(1).EQ.0) GO TO 1015
      L=L+1
      DWISE(KBLOCK,L)=WISE(IBLOCK,2)/PWISE(1)
 1015 IF(KI(2).EQ.0) GO TO 1016
      L=L+1
      DWISE(KBLOCK,L)=1.0
C.....LINAC ALTERATION FINISHING POINT
 1016 WISE(IBLOCK,2)=WISE(IBLOCK,2)+BACKGD
      IF(WISE(IBLOCK,1).LT.YMASK) GO TO 1025
      DY=WISE(IBLOCK,1)-WISE(IBLOCK,2)
      ABSDY=ABS (DY)
      RZNUM=RZNUM+ABSDY
      GO TO (1017,1018,1019,1020),IW
 1017 WEIGHT=1.0
      GO TO 1021
 1018 WEIGHT=1.0/WISE(IBLOCK,1)
      GO TO 1021
 1019 WEIGHT=1.0/(WISE(IBLOCK,1)**2)
      GO TO 1021
 1020 WEIGHT=ABSDY
      GO TO 1021
C     WEIGHT IS THE SQUARE OF ORFLS WEIGHT
C     WDY HAS BEEN REMOVED
 1021 SIG=SIG+WEIGHT*DY**2
      JK=1
      DO 1024 J=1,NV
      TEMP=DWISE(KBLOCK,J)*WEIGHT
      IF(TEMP.NE.0.0)  GO TO 1022
C     BY-PASS IF DERIVATIVE IS ZERO
      JK=JK+NV+1-J
      GO TO 1024
 1022 DO 1023 K=J,NV
      AM(JK)=AM(JK)+TEMP*DWISE(KBLOCK,K)
      JK=JK+1
 1023 CONTINUE
      V(J)=V(J)+TEMP*DY
 1024 JK=JK+LUM
 1025 CONTINUE
C     END LOOP TO STORE MATRIX AND VECTOR
C     TRANSFER BLOCK
      DO 1026 L = 1,KBL
      N=L+KBL
      DO 1026 M = 1,NV
      DWISE(L,M)=DWISE(N,M)
 1026 DWISE(N,M)=0.0
      NBLOCK=NBLOCK+KBL
      IF(NBLOCK.GT.NMAX) GO TO 1028
      BLOCK=NBLOCK*STEP
      IF(IQUERY.LT.NBRAGG) GO TO 1009
      IF(IQUERY.GT.NBRAGG) GO TO 1027
      IQUERY=IQUERY+1
      GO TO 1014
 1027 CONTINUE
C     END LOOP THROUGH OBSERVATIONS
 1028 IF(IC.EQ.NOUT) WRITE(IOO,44) (WISE(I,1),WISE(I,2),I=MINPT,MAXPT)
      IF(IC.EQ.-NOUT) CALL PLOTLP
C     COMPUTE AND PUT OUT AGREEMENT FACTORS
      SQSIG=SQRT (ABS (SIG/FLOAT (NO-NV+LUM)))
      MEANWD=MEANWD/NBRAGG
      NEWNO=NO/MEANWD
      SQSNY=SQRT (ABS (SIG/FLOAT (NEWNO-NV+LUM)))
      WRITE(IOO,57) MEANWD,NEWNO,SQSNY
      SQSNY=SQSNY/SQSIG
      WRITE(IOO,46) IC,SIG,SQSIG
      RZ=RZNUM/YSUM
      WRITE(IOO,41)
      WRITE(IOO,42) RZNUM,YSUM,RZ
C     START LOOP TO TEST FOR ZERO DIAGONAL ELEMENT
      ISING=0
      II=1
      IF(LUM.NE.0) CALL BORDER
      IID=NVL
      DO 1030 I=1,NVL
      IF(AM(II).NE.0.0)  GO TO 1029
      ISING=1
      WRITE(IOO,48)  I
 1029 II=II+IID
      IID=IID-1
 1030 CONTINUE
C     END LOOP TO TEST FOR ZERO DIAGONAL ELEMENT
C     TERMINATE JOB IF ZERO DIAGONAL ELEMENT WAS FOUND
      IF(ISING.NE.0)  CALL FINISH(11)
C     ENTER SUBROUTINE TO REPLACE MATRIX WITH INVERSE
      CALL MATINV(NVL,ISING)
      IF(ISING.NE.0) CALL FINISH(10)
C     START LOOP FOR MATRIX VECTOR MULTIPLICATION FOR PARAMETER CHANGES
      DO 1033 I=1,NV
      PDI=0.0
      IJ=I
      IJD=NVL-1
      DO 1032 J=1,NVL
      PDI=PDI+AM(IJ)*V(J)
      IF(J.GE.I)  GO TO 1031
      IJ=IJ+IJD
      IJD=IJD-1
      GO TO 1032
 1031 IJ=IJ+1
 1032 CONTINUE
      DC(I)=PDI
      SIG=SIG-PDI*V(I)
 1033 CONTINUE
C     END LOOP FOR MATRIX VECTOR MULTIPLICATION
C     RECOMPUTE AGREEMENT FACTOR USING MODIFIED SIG
      SQSIG=SQRT (ABS (SIG/FLOAT (NO-NV+LUM)))
C     PUT OUT CAPTION FOR LIST OF CORRECTED PARAMETERS
      WRITE(IOO,49) IC
C     START LOOP TO CORRECT AND PUT OUT PARAMETERS
      J=1
      II=1
      IJ=NVL
      WRITE(IOF,50) NV
      DO 1034 I=1,NP
      IF(KI(I).EQ.0) GO TO 1034
      POLD=P(I)
      PDI=DC(J)
      P(I)=POLD+PDI
      DC(J)=AM(II)
      II=II+IJ
      IJ=IJ-1
      SIGP=SQRT(DC(J))*SQSIG
      SIGNEW=SIGP*SQSNY
      SIGRAT=PDI/SIGP
      WRITE(IOO,50) I,POLD,PDI,P(I),SIGP,SIGRAT,SIGNEW
      WRITE(IOF,50) I,POLD,PDI,P(I),SIGP,SIGRAT,SIGNEW
      J=J+1
 1034 CONTINUE
      IF(LATTIS.LT.8) GO TO 1035
      P(11)=P(10)
      IF(LATTIS.LT.12) GO TO 1035
      IF(LATTIS.GT.19) GO TO 1035
      P(12)=P(10)
      P(14)=P(13)
      P(15)=P(13)
C     END LOOP TO CORRECT AND PUT OUT PARAMETERS
C     PUT OUT ESTIMATED AGREEMENT FACTORS
 1035 WRITE(IOO,47) IC,SIG,SQSIG
C     ENTER SUBROUTINE TO TEST PARAMETERS
      ISTOP=0
      CALL TEST
      IF(ISTOP.NE.0) CALL FINISH(9)
 1036 CONTINUE
C     END LOOP THROUGH NC CYCLES AND FINAL CALC OF Y
C     TERMINATE JOB
      IF(NC.EQ.0)  GO TO 1041
C     CALCULATE AND PUT OUT CORRELATION MATRIX
      WRITE(IOO,51)
      DO 1037 I=1,NV
      DC(I)=1.0/SQRT(DC(I))
 1037 CONTINUE
      IJ=1
      DO 1040 I=1,NV
      IF(I.NE.1)  V(I-1)=10000.0
      DO 1038 J=I,NV
      V(J)=0.0
 1038 CONTINUE
      DO 1039 J=I,NV
      V(J)=AM(IJ)*DC(I)*DC(J)
      IJ=IJ+1
 1039 CONTINUE
      K=1+10*((I-1)/10)
      WRITE(IOO,52) I,(V(J),J=K,NV)
 1040 IJ=IJ+LUM
 1041 CALL FINISH(12)
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE CONSTR(INPUT,I)
 1984 FORMAT(' BASIC IDMOL CONSTRAINT PROCEDURE')
      COMMON IOI,IOO,AM(861),V(41),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB,LATTIS,MINRH,MAXRH,NBRAGG,INDEX(5,1500),STEP,
     2 WISE(3100,2),MINPT,MAXPT,PTMIN,PTMAX,RMAX,ISCREW(3),IGLIDE(3),
     3 NREFS,IC,NC,NV,NVM,NP,NM,ITF,ISTOP,KEY,NPC,NA,LVAR(41),TWOPIS,
     4 BETA(6,40),P(400),KI(400),XYZ(3,40),SF(40),AI(40),TO,PWISE(9),
     5 BH,DADAI(40),DADX(3,40),DADB(6,40),H(3),NF,ICENT,NS,NREFL,RHO,
     6 NOUT,NPM,TWOPI,YC,SINTHL,DC(400),FX(32,10),TS(3,24),IS(2,3,24),
     7 NO,NATOM,NNA,NNB,NNC,NND,NNE,NNF,NNG,NNH,LUM,NVL,YSUM,YMASK,
     8 YAV,NFRAC,NMAX,NEXT,IOF,IW,INHKL
C  GENERALISED IDMOL  *****************DIMENSIONS***********************
C     DIMENSION XCGS( X),EULS( X),ER( X,3),EQ( X,9)     X=3 OR 6 OR 9
      COMMON /CON/   XCGS(21),EULS(21),ER(21,3),EQ(21,9)
C     DIMENSION XYZS(3, Y),DXYZS(3, Y)     Y IS MAXIMUM 'NNA'
      COMMON /CON/   XYZS(3,26),DXYZS(3,26)
      COMMON /CON/   EE(3,3),EA(3,3),EB(3,3),EC(3,3),MOLSS,MOLSP
      COMMON /CON/   EU(3,3),EV(3,3),EW(3,3),EX(3,3),EY(3,3),EZ(3,3)
C  GENERALISED IDMOL  *************END DIMENSIONS***********************
      IF(INPUT.EQ.1) GO TO 1
      WRITE(IOO,1984)
      IF(INPUT.EQ.2) RETURN
C  GENERALISED IDMOL  *****************SETUP****************************
C  MOLS=NO. OF MOLECULAR UNITS
C     'NNH' MUST NOT EXCEED 'X'/3, 'X' IN DIMENSIONS
      MOLSS=NNH*3
      MOLSP=MOLSS-2
C  NNA=NO. IN UNIT, NNB=NO. NOT IN UNITS
      DO 7007  L = 1,MOLSS
      EULS(L)=P(I)
      XCGS(L)=P(I+MOLSS)
 7007 I=I+1
      I=I+MOLSS
      DO 9203 K = 1,NNC
      DO 9203 J=1,3
      XYZS(J,K)=P(I)
 9203 I=I+1
      IF(NND.EQ.0) GO TO 9205
      DO 9204 K=1,NND
      II=K+NNC
      XYZS(1,II)= XYZS(1,K)
      XYZS(2,II)=-XYZS(2,K)
      XYZS(3,II)= XYZS(3,K)
 9204 CONTINUE
 9205 CONTINUE
      II=-NNA
      DO 7008 LP=1,MOLSP,3
      II=II+NNA
      SPHI=SIN(EULS(LP  ))
      STHE=SIN(EULS(LP+1))
      SPSI=SIN(EULS(LP+2))
      CPHI=COS(EULS(LP  ))
      CTHE=COS(EULS(LP+1))
      CPSI=COS(EULS(LP+2))
      DO 9200 K = 1,3
      DO 9200 J = 1,3
      EU(K,J)=0.0
      EV(K,J)=0.0
      EW(K,J)=0.0
      EX(K,J)=0.0
      EY(K,J)=0.0
 9200 EZ(K,J)=0.0
      EX(3,3)=1.0
      EY(1,1)=1.0
      EZ(3,3)=1.0
      EX(1,1)= CPHI
      EX(2,2)= CPHI
      EX(1,2)= SPHI
      EX(2,1)=-SPHI
      EY(3,3)= CTHE
      EY(2,2)= CTHE
      EY(2,3)= STHE
      EY(3,2)=-STHE
      EZ(1,1)= CPSI
      EZ(2,2)= CPSI
      EZ(1,2)= SPSI
      EZ(2,1)=-SPSI
      EU(1,1)=-SPHI
      EU(2,2)=-SPHI
      EU(1,2)= CPHI
      EU(2,1)=-CPHI
      EV(3,3)=-STHE
      EV(2,2)=-STHE
      EV(2,3)= CTHE
      EV(3,2)=-CTHE
      EW(1,1)=-SPSI
      EW(2,2)=-SPSI
      EW(1,2)= CPSI
      EW(2,1)=-CPSI
      CALL PROD(EA,EV,EX,3)
      CALL PROD(EB,EZ,EA,3)
      CALL PROD(EV,EZ,EY,3)
      CALL PROD(EE,EV,EX,3)
      CALL PROD(EA,EV,EU,3)
      CALL PROD(EV,EY,EX,3)
      CALL PROD(EC,EW,EV,3)
      CALL PROD(EU,ZT,EE,3)
    9 FORMAT(3(/3F12.8,F15.8,2F12.8),17H  TRANSFORMATIONS/)
      WRITE(IOO,9) ((EE(J,K),K=1,3),(EU(J,K),K=1,3),J=1,3)
      DO 9201 K = 1,NNA
      KB=NNB+K+II
      DO 9201 L = 1,3
 9201 XYZ(L,KB)=XCGS(L+LP-1)+
     1      EE(L,1)*XYZS(1,K)+EE(L,2)*XYZS(2,K)+EE(L,3)*XYZS(3,K)
      DO 7009 K = 1,3
      DO 7009 L = 1,3
      KL=K+3*(L-1)
      LPM=LP+K-1
      ER(LPM,L)=EE(K,L)
      EQ(LP  ,KL)=EA(K,L)
      EQ(LP+1,KL)=EB(K,L)
 7009 EQ(LP+2,KL)=EC(K,L)
 7008 CONTINUE
      GO TO 13
C     JUMP HERE TO AVOID SETUP
   11 I=I+NNC*3+NNH*6
   13 CONTINUE
C  GENERALISED IDMOL  *************END SETUP****************************
      GO TO 2
    1 CONTINUE
C  GENERALISED IDMOL  *****************DIFFERENTIALS********************
C     DIFFERENTIALS TO ZERO AT OUTSET
      DO 7005 L = 1,MOLSS
      EULS(L)=0.0
 7005 XCGS(L)=0.0
      DO 7001 K = 1,NNA
C     RUN OVER ATOMS IN THE BASIC UNIT, LOCAL INDEX
      DO 7002 J = 1,3
C     FOR X,Y,Z OR FOR EULER ANGLES
      OH=0.0
      N=K+NNB
C     RELATES ATOMS IN DIFFERENT UNITS, GLOBAL INDEX
      DO 7003 LP=1,MOLSP,3
C     CYCLE OVER THE DIFFERENT UNITS
      LQ=LP+J-1
C     COMBINES LP AND J LOOPS
      DO 7004 M = 1,3
      L=LP+M-1
C     COMBINES LP AND M LOOPS
      OH=OH+ER(L,J)*DADX(M,N)
C     SUM FOR THE ATOM DIFFERENTIALS
      EULS(LQ)=EULS(LQ)+DADX(M,N)*(  EQ(LQ,M  )*XYZS(1,K)
     2                              +EQ(LQ,M+3)*XYZS(2,K)
     3                              +EQ(LQ,M+6)*XYZS(3,K)  )
C     SUM FOR THE EULER ANGLE DIFFERENTIALS
 7004 CONTINUE
      XCGS(LQ)=XCGS(LQ)+DADX(J,N)
C     SUM FOR MOLECULE CENTRE DIFFERENTIALS
 7003 N=N+NNA
      DXYZS(J,K)=OH
 7002 CONTINUE
 7001 CONTINUE
      DO 7006 L = 1,MOLSS
      DC(I)=EULS(L)
      DC(I+MOLSS)=XCGS(L)
 7006 I=I+1
      I=I+MOLSS
      IF(NND.EQ.0) GO TO 9219
      DO 9215 K=1,NND
      II=K+NNC
      DC(I  )=DXYZS(1,K)+DXYZS(1,II)
      DC(I+1)=DXYZS(2,K)-DXYZS(2,II)
      DC(I+2)=DXYZS(3,K)+DXYZS(3,II)
C     COL 25 IS + FOR M SYMMETRY, - FOR 2 SYMMETRY
 9215 I=I+3
 9219 NNDP=NND+1
      DO 9216 K=NNDP,NNC
      DO 9216 J=1,3
      DC(I)=DXYZS(J,K)
 9216 I=I+1
      GO TO 14
C     JUMP HERE TO AVOID DIFFERENTIALS
   12 I=I+NNC*3+NNH*6
   14 CONTINUE
C  GENERALISED IDMOL  *************END DIFFERENTIALS********************
    2 RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE PRELIM
      COMMON IOI,IOO,AM(861),V(41),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB,LATTIS,MINRH,MAXRH,NBRAGG,INDEX(5,1500),STEP,
     2 WISE(3100,2),MINPT,MAXPT,PTMIN,PTMAX,RMAX,ISCREW(3),IGLIDE(3),
     3 NREFS,IC,NC,NV,NVM,NP,NM,ITF,ISTOP,KEY,NPC,NA,LVAR(41),TWOPIS,
     4 BETA(6,40),P(400),KI(400),XYZ(3,40),SF(40),AI(40),TO,PWISE(9),
     5 BH,DADAI(40),DADX(3,40),DADB(6,40),H(3),NF,ICENT,NS,NREFL,RHO,
     6 NOUT,NPM,TWOPI,YC,SINTHL,DC(400),FX(32,10),TS(3,24),IS(2,3,24),
     7 NO,NATOM,NNA,NNB,NNC,NND,NNE,NNF,NNG,NNH,LUM,NVL,YSUM,YMASK,
     8 YAV,NFRAC,NMAX,NEXT,IOF,IW,INHKL
      DIMENSION TEM(3,3)
   71 FORMAT(6F12.5)
   72 FORMAT(24I3)
   73 FORMAT (8F9.3)
   74 FORMAT(3(F11.6,2I2))
   75 FORMAT(//I10,' = NUMBER OF REFLECTIONS'/I10,' = NUMBER OF PARAMETE
     1RS READ')
   76 FORMAT(10F8.0)
   77 FORMAT(7I4)
   78 FORMAT(I10,28H = NUMBER OF OBSERVED POINTS/F10.0,15H = TOTAL COUNT
     1S/F10.3,16H = AVERAGE COUNT/)
   79 FORMAT(' REFLECTIONS INPUT'/10(I4,3I3))
   80 FORMAT(3I4,I6)
      READ(IOI,72)  NC,NP,NPC,ITF,NF,NA,ICENT,NS,NFRAC,NOUT,IW,
     1              NNA,NNB,NNC,NND,NNE,NNF,NNG,NNH,LUM,IRHOMB
      WRITE(IOO,72) NC,NP,NPC,ITF,NF,NA,ICENT,NS,NFRAC,NOUT,IW,
     1              NNA,NNB,NNC,NND,NNE,NNF,NNG,NNH,LUM,IRHOMB
      IF(NF.EQ.0)  GO TO 101
      READ(IOI,73)  ((FX(I,J),I=1,32),J=1,NF)
      WRITE(IOO,73) ((FX(I,J),I=1,32),J=1,NF)
  101 READ(IOI,74)  ((TS(I,J),(IS(K,I,J),K=1,2),I=1,3),J=1,NS)
      WRITE(IOO,74) ((TS(I,J),(IS(K,I,J),K=1,2),I=1,3),J=1,NS)
C     SET UP TRANSFORMATION TO ORTHOGONAL AXES
      READ(IOI,71)  (RCP(I),I=1,6)
      CALL CELLAZ
      READ(IOI,71)  (PWISE(I),I=1,9)
      WRITE(IOO,71) (PWISE(I),I=1,9)
      READ(IOI,71)  TO,BH
      WRITE(IOO,71) TO,BH
      DO 111 I = 1,NA
      READ(IOI,71)  SF(I),AI(I),(XYZ(K,I),K=1,3),BETA(1,I)
      WRITE(IOO,71) SF(I),AI(I),(XYZ(K,I),K=1,3),BETA(1,I)
      IF(NFRAC.NE.1) GO TO 104
      DO 102 J = 1,3
  102 TEM(1,J)=AT(J,1)*XYZ(1,I)+AT(J,2)*XYZ(2,I)+AT(J,3)*XYZ(3,I)
      DO 103 J = 1,3
  103 XYZ(J,I)=TEM(1,J)
  104 IF(I.LE.ITF) GO TO 111
      READ(IOI,71)  (BETA(J,I),J=1,6)
      WRITE(IOO,71) (BETA(J,I),J=1,6)
      IF(NFRAC.NE.1) GO TO 109
      TEM(1,1)=BETA(1,I)
      TEM(2,2)=BETA(2,I)
      TEM(3,3)=BETA(3,I)
      TEM(2,3)=BETA(4,I)
      TEM(3,2)=BETA(4,I)
      TEM(3,1)=BETA(5,I)
      TEM(1,3)=BETA(5,I)
      TEM(1,2)=BETA(6,I)
      TEM(2,1)=BETA(6,I)
      DO 108 J = 1,3
      DO 108 K = 1,3
      TEMP=0.0
      DO 107 L = 1,3
      DO 107 M = 1,3
  107 TEMP=TEMP+AT(J,L)*AT(K,M)*TEM(L,M)
      N=J
      IF(J.NE.K) N=9-J-K
  108 BETA(N,I)=TEMP
      GO TO 111
  109 DO 110 J = 1,6
  110 BETA(J,I)=BETA(J,I)*TWOPIS
  111 CONTINUE
      N=NPC+17+NA*11-ITF*5
C     17 ARE THE 15 + BH AND OVERALL B
      IF(NP.NE.N) CALL FINISH(3)
      I=NP-NPC+1
      IF(NPC.NE.0) READ(IOI,71)  (P(J),J=I,NP)
      IF(NPC.NE.0) WRITE(IOO,71) (P(J),J=I,NP)
      IF(NP.GT.NPM) CALL FINISH(4)
      IF(NA.GT.NATOM) CALL FINISH(5)
C     INPUT OF REFLECTIONS AND SCAN
      READ(IOI,77) LATTIS,(ISCREW(I),I=1,3),(IGLIDE(I),I=1,3)
      IF(MINPT.EQ.0) CALL FINISH(13)
      INHKL=0
      IF(NBRAGG.NE.0) GO TO 112
      CALL REFGEN
      GO TO 113
  112 READ(IOI,80)  ((INDEX(I,J),I=1,4),J=1,NBRAGG)
      WRITE(IOO,79) ((INDEX(I,J),I=1,4),J=1,NBRAGG)
      INHKL=1
  113 WRITE(IOO,75) NBRAGG,NP
      DO 114 I = 1,NMAX
  114 WISE(I,1)=YMASK-1.
      READ(IOI,76) (WISE(I,1),I=MINPT,MAXPT)
      YSUM=0.0
      NO=0
      DO 115 I=MINPT,MAXPT
      IF(WISE(I,1).LT.YMASK) GO TO 115
      NO=NO+1
      YSUM=YSUM+WISE(I,1)
  115 CONTINUE
      YAV=YSUM/NO
      WRITE(IOO,78) NO,YSUM,YAV
      IF(MAXPT.GT.NMAX) CALL FINISH(7)
      IC=0
      CALL TEST
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE REFGEN
      COMMON IOI,IOO,AM(861),V(41),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB,LATTIS,MINRH,MAXRH,NBRAGG,INDEX(5,1500),STEP,
     2 WISE(3100,2),MINPT,MAXPT,PTMIN,PTMAX,RMAX,ISCREW(3),IGLIDE(3),
     3 NREFS
C SEE PAGE 32 OF INTERNATIONAL TABLES, VOLUME 1
C  LATT = 0  1  2  3  4  5            PROGRAM DEDUCES 'LATT'
C         P  R  C  F  I  H  LATTICE
C VALUE FOR R AND H LATTICES DEPENDS ON POINT GROUP
C LAUE GROUP    LATTIS NUMBERS
C               P    C    F    I    LATTICES
C -1            1
C 2/M  Y UNIQUE 2    3
C MMM           4    5    6    7
C 4/M           8              9
C 4/MMM         10             11
C M3            12        13   14
C M3M           15        16   17
C  RHOMB 18 R-3  R3  R-3C
C  RHOMB 19 R32  R3M  R3C  R-3M
C -3    HEX 20
C -3M   HEX 21
C 6/M   HEX 22
C 6/MMM HEX 23
C SCREW AND GLIDE ABSENCES REMOVED, ORTHORHOMBIC AND MONOCLINIC
C BY SETTING THE MULTIPLICITY TO ZERO
C ISCREW(I) =0 : NO SCREW AXIS
C ISCREW(I) =I : 2 FOLD SCREW AXIS ALONG AXIS I
C IGLIDE(I) =0 : NO GLIDE PLANE PERPENDICULAR TO AXIS I
C  "        =1 : A GLIDE PERPENDICULAR TO AXIS I
C  "        =2 : B GLIDE PERPENDICULAR TO AXIS I
C  "        =3 : C GLIDE PERPENDICULAR TO AXIS I
C  "        =4 : N GLIDE PERPENDICULAR TO AXIS I
 3008 FORMAT(//' REFLECTIONS TO BE GENERATED'/' LATTICE TYPE',I4)
 3011 FORMAT(' 2-FOLD SCREW AXIS ALONG AXIS',I2)
 3012 FORMAT(' GLIDE PLANE PERPENDICULAR TO AXIS',I2,' ALONG AXIS',I2)
 3013 FORMAT(' N GLIDE PERPENDICULAR TO AXIS',I2)
 3014 FORMAT(' MAXIMUM HKL = ',I4)
      DO 1008 I=1,NREFS
      DO 1008 J=1,5
 1008 INDEX(J,I)=0
      WRITE(IOO,3008) LATTIS
      RCPMAX=0.0
      DO 1010 I = 1,3
      IF(ISCREW(I).EQ.I) WRITE(IOO,3011) I
      DO 1009 J = 1,3
      IF(IGLIDE(I).EQ.J) WRITE(IOO,3012) I,J
 1009 CONTINUE
      IF(IGLIDE(I).EQ.4) WRITE(IOO,3013) I
      IF(RCP(I).GT.RCPMAX) RCPMAX=RCP(I)
 1010 CONTINUE
      LIMHKL=2.0*(RCPMAX*RMAX+1.5)
      WRITE(IOO,3014) LIMHKL
      LATT=0
      LAUE=0
      NBRAGG=1
      IF(LATTIS.GT.11) GO TO 1022
      IF(LATTIS.LT.8 ) LATT=LATTIS-3
      IF(LATTIS.LT.5 ) LATT=0
      IF(LATTIS.EQ.3 ) LATT=2
      IF(LATTIS.EQ.9 ) LATT=4
      IF(LATTIS.EQ.11) LATT=4
      DO 1016 I = 1,LIMHKL
      CALL SETHKL(LATT,I,0,0,2)
      CALL SETHKL(LATT,0,I,0,2)
      CALL SETHKL(LATT,0,0,I,2)
      DO 1016 J = 1,LIMHKL
      IF(LATTIS.EQ.1) GO TO 1011
      CALL SETHKL(LATT,I,J,0,4)
      CALL SETHKL(LATT,0,I,J,4)
      IF(LATTIS.GT.3) GO TO 1012
 1011 CALL SETHKL(LATT,I,0, J,2)
      CALL SETHKL(LATT,I,0,-J,2)
      IF(LATTIS.NE.1) GO TO 1013
      CALL SETHKL(LATT,I, J, 0,2)
      CALL SETHKL(LATT,I,-J, 0,2)
      CALL SETHKL(LATT,0, I, J,2)
      CALL SETHKL(LATT,0, I,-J,2)
      GO TO 1013
 1012 CALL SETHKL(LATT,I,0,J,4)
 1013 DO 1016 K = 1,LIMHKL
      IF(LATTIS.LT.4) GO TO 1014
      CALL SETHKL(LATT,I,J,K,8)
      GO TO 1016
 1014 IF(LATTIS.EQ.1) GO TO 1015
      CALL SETHKL(LATT,I,J, K,4)
      CALL SETHKL(LATT,I,J,-K,4)
      GO TO 1016
 1015 CALL SETHKL(LATT,I, J, K,2)
      CALL SETHKL(LATT,I, J,-K,2)
      CALL SETHKL(LATT,I,-J, K,2)
      CALL SETHKL(LATT,I,-J,-K,2)
 1016 CONTINUE
      IF(LATTIS.LT.8) GO TO 1033
C ORTHORHOMBIC, MONOCLINIC AND TRICLINIC DONE ........
      IF(LATTIS.LT.10) GO TO 1019
      DO 1017 I = 2,NBRAGG
      IF(INDEX(4,I).NE.8) GO TO 1017
      IF(INDEX(1,I).GT.INDEX(2,I)) INDEX(4,I)=16
      IF(INDEX(1,I).LT.INDEX(2,I)) INDEX(4,I)=0
 1017 CONTINUE
      DO 1018 I = 2,NBRAGG
      IF(INDEX(4,I).NE.4) GO TO 1018
      IF(INDEX(3,I).NE.0) GO TO 1018
      IF(INDEX(1,I).GT.INDEX(2,I)) INDEX(4,I)=8
      IF(INDEX(1,I).LT.INDEX(2,I)) INDEX(4,I)=0
 1018 CONTINUE
C 4/MMM TETRAG. PART DONE ........
C 4/M AND 4/MMM TO FINISH ........
 1019 DO 1020 I = 2,NBRAGG
      IF(INDEX(4,I).NE.2) GO TO 1020
      IF(INDEX(1,I).NE.0) INDEX(4,I)=4
      IF(INDEX(2,I).NE.0) INDEX(4,I)=0
 1020 CONTINUE
      DO 1021 I = 2,NBRAGG
      IF(INDEX(4,I).NE.4) GO TO 1021
      IF(INDEX(3,I).EQ.0) GO TO 1021
      INDEX(4,I)=8
      IF(INDEX(2,I).EQ.0) INDEX(4,I)=0
 1021 CONTINUE
      GO TO 1033
C TETRAGONAL DONE ........
 1022 IF(LATTIS.GT.17) GO TO 1024
      IF(LATTIS.NE.12) LATT=LATTIS-10
      IF(LATTIS.GT.15) LATT=LATTIS-13
      IF(LATTIS.LT.15) LAUE=1
      DO 1023 I = 1,LIMHKL
      CALL SETHKL(LATT,I,0,0,6)
      CALL SETHKL(LATT,I,I,0,12)
      CALL SETHKL(LATT,I,I,I,8)
      L=I+1
      DO 1023 J = L,LIMHKL
      IF(LAUE.EQ.0) CALL SETHKL(LATT,I,J,0,24)
      IF(LAUE.EQ.1) CALL SETHKL(LATT,I,J,0,12)
      IF(LAUE.EQ.1) CALL SETHKL(LATT,J,I,0,12)
      CALL SETHKL(LATT,I,I,J,24)
      CALL SETHKL(LATT,J,J,I,24)
      M=J+1
      DO 1023 K = M,LIMHKL
      IF(LAUE.EQ.0) CALL SETHKL(LATT,I,J,K,48)
      IF(LAUE.EQ.1) CALL SETHKL(LATT,I,J,K,24)
      IF(LAUE.EQ.1) CALL SETHKL(LATT,K,J,I,24)
 1023 CONTINUE
      GO TO 1033
C CUBIC DONE ........
 1024 IF(LATTIS.GT.19) GO TO 1028
      N=12
      IF(LATTIS.EQ.18) N=6
      DO 1027 I = 1,LIMHKL
      CALL SETHKL(0,I,0,0,6)
      CALL SETHKL(0,I,I,0,6)
      CALL SETHKL(0,I,-I,0,6)
      CALL SETHKL(0,I,I,I,2)
      CALL SETHKL(0,I,I,-I,6)
      L=I+1
      DO 1027 J = L,LIMHKL
      CALL SETHKL(0,J,I,I,6)
      CALL SETHKL(0,-J,I,I,6)
      CALL SETHKL(0,J,J,I,6)
      CALL SETHKL(0,J,J,-I,6)
      IF(N.EQ.12) GO TO 1025
      CALL SETHKL(0,J,I,0,6)
      CALL SETHKL(0,-J,I,0,6)
      CALL SETHKL(0,-I,I,J,6)
      CALL SETHKL(0,-J,J,I,6)
 1025 CALL SETHKL(0,I,J,0,N)
      CALL SETHKL(0,I,-J,0,N)
      CALL SETHKL(0,I,-I,J,N)
      CALL SETHKL(0,J,-J,I,N)
      M=J+1
      DO 1027 K = M,LIMHKL
      IF(N.EQ.12) GO TO 1026
      CALL SETHKL(0,J,K,I,6)
      CALL SETHKL(0,J,K,-I,6)
      CALL SETHKL(0,-J,K,I,6)
      CALL SETHKL(0,J,-K,I,6)
 1026 CALL SETHKL(0,K,J,I,N)
      CALL SETHKL(0,K,J,-I,N)
      CALL SETHKL(0,K,-J,I,N)
 1027 CALL SETHKL(0,-K,J,I,N)
      GO TO 1033
C RHOMBOHEDRAL DONE ........
 1028 N=0
      IF(LATTIS.LT.22) N=5
      M=0
      IF(LATTIS.EQ.20) M=5
      IF(LATTIS.EQ.20) LAUE=1
      IF(LATTIS.EQ.22) LAUE=1
      DO 1032 I = 1,LIMHKL
      L=I+1
      CALL SETHKL(0,I,0,0,6)
      CALL SETHKL(0,I,I,0,6)
      CALL SETHKL(0,0,0,I,2)
      DO 1030 K = 1,LIMHKL
      CALL SETHKL(M,I,I,K,12)
      CALL SETHKL(N,I,0,K,12)
      DO 1030 J = L,LIMHKL
      IF(LAUE.EQ.0) GO TO 1029
      CALL SETHKL(N,I,J,K,12)
      CALL SETHKL(N,J,I,K,12)
      GO TO 1030
 1029 CALL SETHKL(N,I,J,K,24)
 1030 CONTINUE
      DO 1032 J = 1,LIMHKL
      IF(LAUE.EQ.0) GO TO 1031
      CALL SETHKL(0,I,J,0,6)
      CALL SETHKL(0,J,I,0,6)
      GO TO 1032
 1031 CALL SETHKL(0,I,J,0,12)
 1032 CONTINUE
C HEXAGONAL LATTICES DONE ........
C REMOVE SCREW AND GLIDE ABSENCES ........
 1033 IF(NBRAGG.EQ.NREFS) NREFS=-1
      DO 1037 L = 1,3
      M=L+1
      IF(M.EQ.4) M=1
      N=6-L-M
      J=ISCREW(L)
      K=IGLIDE(L)
      IF(J.EQ.0) GO TO 1035
      DO 1034 I = 2,NBRAGG
      IF(INDEX(M,I).NE.0.OR.INDEX(N,I).NE.0) GO TO 1034
      IPAR=INDEX(L,I)
      IF(MOD(IPAR,2).NE.0) INDEX(4,I)=0
 1034 CONTINUE
 1035 IF(K.EQ.0) GO TO 1037
      DO 1036 I = 2,NBRAGG
      IF(INDEX(L,I).NE.0) GO TO 1036
      IPAR=INDEX(K,I)
      IF(K.EQ.4) IPAR=INDEX(M,I)+INDEX(N,I)
      IF(MOD(IPAR,2).NE.0) INDEX(4,I)=0
 1036 CONTINUE
 1037 CONTINUE
C REMOVE THE ZERO MULTIPLICITIES ........
      N=1
      DO 1039 I = 2,NBRAGG
      IF(INDEX(4,I).EQ.0) GO TO 1039
      N=N+1
      DO 1038 K = 1,5
 1038 INDEX(K,N)=INDEX(K,I)
 1039 CONTINUE
C SORT ON INCREASING RHO ........
      NBRAGG=N-1
      DO 1041 I = 1,NBRAGG
      L=I+1
      IRHO=MAXRH
      DO 1040 J = L,N
      IF(INDEX(5,J).GT.IRHO) GO TO 1040
      IRHO=INDEX(5,J)
      M=J
 1040 CONTINUE
      DO 1041 K = 1,5
      INDEX(K,I)=INDEX(K,M)
 1041 INDEX(K,M)=INDEX(K,L)
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE SETHKL(LATT,I,J,K,M)
      COMMON IOI,IOO,AM(861),V(41),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB,LATTIS,MINRH,MAXRH,NBRAGG,INDEX(5,1500),STEP,
     2 WISE(3100,2),MINPT,MAXPT,PTMIN,PTMAX,RMAX,ISCREW(3),IGLIDE(3),
     3 NREFS
      DIMENSION HJ(3),HK(3)
      IF(NBRAGG.EQ.NREFS) RETURN
      IF(LATT.LT.2) GO TO 1043
      IF(LATT.GT.4) GO TO 1043
      L=I+J
      IF(LATT.EQ.4) L=L+K
      IF(MOD(L,2).NE.0) RETURN
      IF(LATT.NE.3) GO TO 1043
      L=I+K
      IF(MOD(L,2).NE.0) RETURN
 1043 HK(1)=I
      HK(2)=J
      HK(3)=K
      HJ(1)=ZT(1,1)*HK(1)+ZT(2,1)*HK(2)+ZT(3,1)*HK(3)
      HJ(2)=ZT(1,2)*HK(1)+ZT(2,2)*HK(2)+ZT(3,2)*HK(3)
      HJ(3)=ZT(1,3)*HK(1)+ZT(2,3)*HK(2)+ZT(3,3)*HK(3)
      IRHO=2500.*(HJ(1)**2+HJ(2)**2+HJ(3)**2)
      IF(IRHO.GT.MAXRH) RETURN
      IF(IRHO.LT.MINRH) RETURN
      NBRAGG=NBRAGG+1
      INDEX(1,NBRAGG)=I
      INDEX(2,NBRAGG)=J
      INDEX(3,NBRAGG)=K
      INDEX(4,NBRAGG)=M
      INDEX(5,NBRAGG)=IRHO
      IF(LATT.NE.5) RETURN
      INDEX(4,NBRAGG)=M/2
      IF(NBRAGG.EQ.NREFS) RETURN
      NBRAGG=NBRAGG+1
      INDEX(1,NBRAGG)=I
      INDEX(2,NBRAGG)=J
      INDEX(3,NBRAGG)=-K
      INDEX(4,NBRAGG)=M/2
      INDEX(5,NBRAGG)=IRHO
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE CALC
      COMMON IOI,IOO,AM(861),V(41),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB,LATTIS,MINRH,MAXRH,NBRAGG,INDEX(5,1500),STEP,
     2 WISE(3100,2),MINPT,MAXPT,PTMIN,PTMAX,RMAX,ISCREW(3),IGLIDE(3),
     3 NREFS,IC,NC,NV,NVM,NP,NM,ITF,ISTOP,KEY,NPC,NA,LVAR(41),TWOPIS,
     4 BETA(6,40),P(400),KI(400),XYZ(3,40),SF(40),AI(40),TO,PWISE(9),
     5 BH,DADAI(40),DADX(3,40),DADB(6,40),H(3),NF,ICENT,NS,NREFL,RHO,
     6 NOUT,NPM,TWOPI,YC,SINTHL,DC(400),FX(32,10),TS(3,24),IS(2,3,24),
     7 NO,NATOM,NNA,NNB,NNC,NND,NNE,NNF,NNG,NNH,LUM,NVL,YSUM,YMASK,
     8 YAV,NFRAC,NMAX,NEXT,IOF,IW,INHKL
      DIMENSION DBDAI(40),DBDX(3,40),DBDB(6,40),HK(3),HJ(3),HHJ(6)
C     COMPUTE REQUIRED SUMS OVER NS SYMMETRY TRANSFORMATIONS
C     FOR EACH OF NA ATOMS
      H(1)=INDEX(1,NREFL)
      H(2)=INDEX(2,NREFL)
      H(3)=INDEX(3,NREFL)
C     CLEAR STORAGE ARRAYS
      DO 3002 I=1,NA
      DADAI(I)=0.0
      DBDAI(I)=0.0
      DO 3001 K=1,3
      DADX(K,I)=0.0
 3001 DBDX(K,I)=0.0
      DO 3002 K=1,6
      DADB(K,I)=0.0
 3002 DBDB(K,I)=0.0
C     START LOOP THROUGH THE NS SYMMETRY TRANSFORMATIONS
      DO 3015 J=1,NS
C     STORE TRANSFORMED INDICES
      TJ=0.0
      DO 3005 I=1,3
      TJ=TJ+H(I)*TS(I,J)
      HK(I)=0.0
      DO 3005 K=1,3
      DO 3004 L=1,2
      M=IS(L,K,J)
      IF(M.NE.I)  GO TO 3003
      HK(I)=HK(I)+H(K)
      GO TO 3004
 3003 M=-M
      IF(M.NE.I)  GO TO 3004
      HK(I)=HK(I)-H(K)
 3004 CONTINUE
 3005 CONTINUE
C     TRANSFORM THE NEW INDICES WITH TRANSPOSED INVERSE MATRIX
      HJ(1)=ZT(1,1)*HK(1)+ZT(2,1)*HK(2)+ZT(3,1)*HK(3)
      HJ(2)=ZT(1,2)*HK(1)+ZT(2,2)*HK(2)+ZT(3,2)*HK(3)
      HJ(3)=ZT(1,3)*HK(1)+ZT(2,3)*HK(2)+ZT(3,3)*HK(3)
C     COMPUTE PRODUCTS OF TRANSFORMED INDICES
      HHJ(1)=HJ(1)**2
      HHJ(2)=HJ(2)**2
      HHJ(3)=HJ(3)**2
      HHJ(4)=2.0*HJ(2)*HJ(3)
      HHJ(5)=2.0*HJ(3)*HJ(1)
      HHJ(6)=2.0*HJ(1)*HJ(2)
C     START LOOP THROUGH THE NA ATOMS
      DO 3014 I=1,NA
C     ENTER USERS SUBROUTINE TO SET INDICES FOR SPECIAL POSITIONS
C     COMPUTE TRIGONOMETRIC TERMS
      HX=TJ
      DO 3006 K=1,3
 3006 HX=HX+HJ(K)*XYZ(K,I)
      HX=HX*TWOPI
      COSIJ=COS(HX)
      SINIJ=SIN(HX)
C     COMPUTE ANISOTROPIC TEMPERATURE FACTOR
      IF(I.LE.ITF) GO TO 3008
      HHB=0.0
      DO 3007 K=1,6
 3007 HHB=HHB+HHJ(K)*BETA(K,I)
      EXPIJ=EXP (-HHB)
      COSIJ=COSIJ*EXPIJ
      SINIJ=SINIJ*EXPIJ
C     ACCUMULATE REQUIRED SUMS
 3008 DADAI(I)=DADAI(I)+COSIJ
      DO 3009 K=1,3
 3009 DADX(K,I)=DADX(K,I)+HJ(K)*SINIJ
      IF(I.LE.ITF) GO TO 3011
      DO 3010 K=1,6
 3010 DADB(K,I)=DADB(K,I)+HHJ(K)*COSIJ
 3011 IF(ICENT.EQ.1)  GO TO 3014
C     BYPASS REMAINING SUMS IF CENTROSYMMETRIC
      DBDAI(I)=DBDAI(I)+SINIJ
      DO 3012 K=1,3
 3012 DBDX(K,I)=DBDX(K,I)+HJ(K)*COSIJ
      IF(I.LE.ITF) GO TO 3014
      DO 3013 K=1,6
 3013 DBDB(K,I)=DBDB(K,I)+HHJ(K)*SINIJ
 3014 CONTINUE
C     END LOOP THROUGH THE NA ATOMS
 3015 CONTINUE
C     END LOOP THROUGH THE NS SYMMETRY TRANSFORMATIONS
C     COMPUTE A AND B AND THEIR DERIVATIVES
      RHO=(HJ(1)**2+HJ(2)**2+HJ(3)**2)*0.25
      INDEX(5,NREFL)=10000.*RHO
      SINTHL=SQRT (RHO)
C     EXTINCTION PART REMOVED
      A=0.0
      B=0.0
      DAH=0.0
      DBH=0.0
      PT=SINTHL/.05+1.0
      IPT=PT
      FRACPT=PT-FLOAT (IPT)
      LAST=0
C     START LOOP THROUGH NA ATOMS
      FI=1.0
      DO 3026 I=1,NA
C     OBTAIN SCATTERING FACTOR
      LET=0
      IF(SF(I).EQ.SF(1)) LET=1
      IF(NF.EQ.0)  GO TO 3016
      L=SF(I)
      IF(L.EQ.LAST)  GO TO 3016
      LAST=L
      FI=FX(IPT,L)+(FX(IPT+1,L)-FX(IPT,L))*FRACPT+FX(32,L)
C     COMPUTE ISOTROPIC TEMPERATURE FACTOR
 3016 IF(I.LE.ITF) GO TO 3017
      TFI=1.0
      GO TO 3018
 3017 TFI=EXP (-BETA(1,I)*RHO)
C     COMPUTE A AND ITS DERIVATIVES
 3018 FTA=FI*TFI*AI(I)
      IF(LET.EQ.1) FTA=FTA*EXP(-RHO*BH)
      FTACOS=FTA*DADAI(I)
      A=A+FTACOS
      DADAI(I)=FI*TFI*DADAI(I)
      PIFI=-TWOPI*FTA
      DO 3019 K=1,3
 3019 DADX(K,I)=PIFI*DADX(K,I)
      IF(LET.EQ.1) DAH=DAH-RHO*FTACOS
      IF(I.LE.ITF)  GO TO 3021
      DO 3020 K=1,6
 3020 DADB(K,I)=-FTA*DADB(K,I)
      GO TO 3022
 3021 DADB(1,I)=-RHO*FTACOS
C     COMPUTE B AND ITS DERIVATIVES UNLESS CENTROSYMMETRIC
 3022 IF(ICENT.EQ.1)  GO TO 3026
      FTASIN=FTA*DBDAI(I)
      B=B+FTASIN
      DBDAI(I)=FI*TFI*DBDAI(I)
      PIFI=-PIFI
      DO 3023 K=1,3
 3023 DBDX(K,I)=PIFI*DBDX(K,I)
      IF(LET.EQ.1) DBH=DBH-RHO*FTASIN
      IF(I.LE.ITF)  GO TO 3025
      DO 3024 K=1,6
 3024 DBDB(K,I)=-FTA*DBDB(K,I)
      GO TO 3026
 3025 DBDB(1,I)=-RHO*FTASIN
 3026 CONTINUE
C     END LOOP THROUGH NA ATOMS
C     COMPUTE SCALE FACTOR TIMES OVERALL TEMPERATURE FACTOR
C     CALCULATE F**2 FOR CENTROSYMMETRIC OR NON-CENTROSYMMETRIC
C     STRUCTURES
      SQTO=PWISE(1)*EXP (-2.0*TO*RHO)
      IF(ICENT.NE.1)  GO TO 3027
      CA=8.0*SQTO*A
      CB=0.0
      YC=CA*A/2.0
      GO TO 3028
 3027 YC=SQTO*(A**2+B**2)
      CA=2.0*SQTO*A
      CB=2.0*SQTO*B
 3028 DFDTO=-2.0*RHO*YC
C     COMPUTE DERIVATIVES OF F**2
      DAH=CA*DAH+CB*DBH
C     START LOOP THROUGH NA ATOMS
      DO 3032 K=1,NA
      DADAI(K)=CA*DADAI(K)+CB*DBDAI(K)
      DO 3029 I=1,3
 3029 DADX(I,K)=CA*DADX(I,K)+CB*DBDX(I,K)
      IF(K.GT.ITF) GO TO 3030
      DADB(1,K)=CA*DADB(1,K)+CB*DBDB(1,K)
      GO TO 3032
 3030 DO 3031 I=1,6
 3031 DADB(I,K)=CA*DADB(I,K)+CB*DBDB(I,K)
 3032 CONTINUE
C     COPY DERIVATIVES INTO ARRAY DC
      I=16
      DC(I)=DFDTO
      DC(I+1)=DAH
      I=I+2
      DO 3036 J=1,NA
      DC(I+1)=DADAI(J)
      I=I+2
      DO 3033 K=1,3
      DC(I)=DADX(K,J)
 3033 I=I+1
      IF(J.GT.ITF) GO TO 3034
      DC(I)=DADB(1,J)
      I=I+1
      GO TO 3036
 3034 DO 3035 K=1,6
      DC(I)=DADB(K,J)
 3035 I=I+1
 3036 CONTINUE
C********** ALTERATION POINT 2 ---- DERIVATIVES, INPUT=1 ***************
      IF(NPC.NE.0) CALL CONSTR(1,I)
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE TEST
      COMMON IOI,IOO,AM(861),V(41),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB,LATTIS,MINRH,MAXRH,NBRAGG,INDEX(5,1500),STEP,
     2 WISE(3100,2),MINPT,MAXPT,PTMIN,PTMAX,RMAX,ISCREW(3),IGLIDE(3),
     3 NREFS,IC,NC,NV,NVM,NP,NM,ITF,ISTOP,KEY,NPC,NA,LVAR(41),TWOPIS,
     4 BETA(6,40),P(400),KI(400),XYZ(3,40),SF(40),AI(40),TO,PWISE(9),
     5 BH,DADAI(40),DADX(3,40),DADB(6,40),H(3),NF,ICENT,NS,NREFL,RHO,
     6 NOUT,NPM,TWOPI,YC,SINTHL,DC(400),FX(32,10),TS(3,24),IS(2,3,24),
     7 NO,NATOM,NNA,NNB,NNC,NND,NNE,NNF,NNG,NNH,LUM,NVL,YSUM,YMASK,
     8 YAV,NFRAC,NMAX,NEXT,IOF,IW,INHKL
      DIMENSION TABLE(14),XYZF(3,40)
   70 FORMAT(72I1)
   72 FORMAT('0TEMPERATURE FACTOR OF ATOM'I3' ISNT POSITIVE-DEFINITE')
   73 FORMAT(2F8.4,2X,3F8.4,2X,3F9.5,2X,6F8.4)
   74 FORMAT(//' PARAMETERS ORTHOGONAL THEN FRACTIONAL AND B(I,J)'/)
   75 FORMAT(/I10,27H = PARAMETERS TO BE VARIED ,(/37X,36I2))
   76 FORMAT(8X,I6,26X,F13.7)
   77 FORMAT(' NUMBER OF LAGRANGE MULTIPLIERS ='I3)
   78 FORMAT(I8,' PARAMETERS AND THEIR INDICES INPUT FROM LAST CYCLE')
      IF(NEXT.NE.0) GO TO 312
      GO TO 303
  301 READ(IOI,76)  NEXT
      WRITE(IOO,78) NEXT
      DO 302 J = 1,NEXT
      READ(IOI,76)  NPAR,PAR
      WRITE(IOO,76) NPAR,PAR
  302 P(NPAR)=PAR
      IC=-1
  303 NEXT=0
      IF(IC.EQ.NC) GO TO 305
      READ(IOI,70) (KI(I),I=1,NP)
      KEY=KI(1)+KI(2)+KI(3)+KI(4)+1
C     SET UP LVAR( ) FOR VARIABLES EXCEPT SCALE AND BACKGROUND
      NV=0
      DO 304 J=1,NP
      IF(KI(J).NE.1) GO TO 304
      NV=NV+1
      LVAR(NV)=J
  304 CONTINUE
      WRITE(IOO,75)  NV,(KI(I),I=1,NP)
      NVL=NV+LUM
      IF(LUM.NE.0) WRITE(IOO,77) LUM
      NM=(NVL*(NVL+1))/2
  305 I=NP-NPC+1
      IF(IC.EQ.0) GO TO 311
      DO 306 J=1,9
      PWISE(J)=P(J)
      IF(J.GT.3) GO TO 306
      RCP(J  )=P(J+9 )
      RCP(J+3)=P(J+12)/DEGRAD
  306 CONTINUE
      TO=P(16)
      BH=P(17)
      I=18
      DO 310 J=1,NA
      SF(J)=P(I)
      AI(J)=P(I+1)
      I=I+2
      DO 307 K=1,3
      XYZ(K,J)=P(I)
  307 I=I+1
      IF(J.GT.ITF) GO TO 308
      BETA(1,J)=P(I)
      I=I+1
      GO TO 310
  308 DO 309 K=1,6
      BETA(K,J)=P(I)
  309 I=I+1
  310 CONTINUE
      CALL PROD(XYZF,ZT,XYZ,NA)
C     END OLD PTOSB
      CALL CELLAZ
      CALL PROD(XYZ,AT,XYZF,NA)
C********** ALTERATION POINT 1 ---- INPUT=0 ****************************
  311 IF(NPC.NE.0) CALL CONSTR(0,I)
      WRITE(IOO,74)
  312 DO 313 J = 1,9
      P(J)=PWISE(J)
      IF(J.GT.3) GO TO 313
      P(J+9 )=RCP(J  )
      P(J+12)=RCP(J+3)*DEGRAD
  313 CONTINUE
      P(16)=TO
      P(17)=BH
      CALL PROD(XYZF,ZT,XYZ,NA)
      I=18
      DO 317 J=1,NA
      TABLE(1)=SF(J)
      TABLE(2)=AI(J)
      P(I)=SF(J)
      P(I+1)=AI(J)
      I=I+2
      DO 314 K=1,3
      P(I)=XYZ(K,J)
      TABLE(K+2)=P(I)
      TABLE(K+5)=XYZF(K,J)
  314 I=I+1
      IF(J.GT.ITF) GO TO 315
      P(I)=BETA(1,J)
      TABLE(9)=P(I)
      I=I+1
      IF(NEXT.EQ.0) WRITE(IOO,73) (TABLE(K),K=1,9)
      GO TO 317
  315 DO 316 K=1,6
      P(I)=BETA(K,J)
      TABLE(K+8)=BETA(K,J)/TWOPIS
  316 I=I+1
      IF(NEXT.EQ.0) WRITE(IOO,73) (TABLE(K),K=1,14)
  317 CONTINUE
      IF(NEXT.NE.0) GO TO 301
      IF(NV.EQ.0) CALL FINISH(1)
      IF(NVL.GT.NVM) CALL FINISH(8)
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE CELLAZ
      COMMON IOI,IOO,AM(861),V(41),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB
      DIMENSION RECIP(6),SN(3),CS(3),RA(3)
   51 FORMAT(//11H NEW CELL =,6F12.6/12H CELL VOLUME ,F15.5/27H RECIPROC
     1AL CELL PARAMETERS,6F10.5/16H TRANSFORMATIONS//(3F12.6,9X,3F12.6))
   52 FORMAT(' SPECIAL RHOMBOHEDRAL SECTION')
      IF(IRHOMB.NE.0) GO TO 1303
      DO 1301 I=1,3
      SN(I)=SIN(RCP(I+3)*DEGRAD)
 1301 CS(I)=COS(RCP(I+3)*DEGRAD)
      RA(1)=CS(2)*CS(3)-CS(1)
      RA(2)=CS(3)*CS(1)-CS(2)
      RA(3)=CS(1)*CS(2)-CS(3)
      BIGD=1.-CS(1)**2-CS(2)**2-CS(3)**2+2.*CS(1)*CS(2)*CS(3)
      AROOT=SQRT(BIGD)
      AT(1,1)=RCP(1)
      AT(1,2)=RCP(2)*CS(3)
      AT(1,3)=RCP(3)*CS(2)
      AT(2,1)=0.
      AT(2,2)=RCP(2)*SN(3)
      AT(2,3)=-RCP(3)*RA(1)/SN(3)
      AT(3,1)=0.
      AT(3,2)=0.
      AT(3,3)=RCP(3)*AROOT/SN(3)
      CELL=AT(1,1)*AT(2,2)*AT(3,3)
      ZT(1,1)=1./RCP(1)
      ZT(1,2)=-CS(3)/(RCP(1)*SN(3))
      ZT(1,3)=RA(2)/(RCP(1)*SN(3)*AROOT)
      ZT(2,1)=0.
      ZT(2,2)=1./(RCP(2)*SN(3))
      ZT(2,3)=RA(1)/(RCP(2)*SN(3)*AROOT)
      ZT(3,1)=0.
      ZT(3,2)=0.
      ZT(3,3)=SN(3)/(RCP(3)*AROOT)
      DO 1302 I = 1,3
      L=I+3
      J=I+1
      IF(J.EQ.4) J=1
      K=6-I-J
      RECIP(I)=SN(I)/(RCP(I)*AROOT)
      RECIP(L)=ARCOS(RA(I)/(SN(J)*SN(K)))/DEGRAD
      FAC=0.5/BIGD
      DCELL(I,I  )=-FAC*SN(I)**2/RCP(I)**3
      DCELL(I,J  )=0.0
      DCELL(I,K  )=0.0
      DCELL(I,I+3)=0.0
      DCELL(I,J+3)=-FAC*RA(J)/(RCP(K)*RCP(I)**2)
      DCELL(I,K+3)=-FAC*RA(K)/(RCP(J)*RCP(I)**2)
      FAC=FAC*SN(I)/BIGD
      DCELL(L,I  )=FAC*(CS(I)*BIGD+SN(I)**2*RA(I))/RCP(I)**2
      DCELL(L,J  )=FAC*            SN(J)**2*RA(I) /RCP(J)**2
      DCELL(L,K  )=FAC*            SN(K)**2*RA(I) /RCP(K)**2
      DCELL(L,I+3)=FAC*(2.*RA(I)*RA(I)+BIGD      )/(RCP(J)*RCP(K))
      DCELL(L,J+3)=FAC*(2.*RA(I)*RA(J)-BIGD*CS(K))/(RCP(K)*RCP(I))
      DCELL(L,K+3)=FAC*(2.*RA(I)*RA(K)-BIGD*CS(J))/(RCP(I)*RCP(J))
 1302 CONTINUE
      GO TO 1308
C     RHOMBOHEDRAL: SPECIAL CASE
 1303 ROOT2=SQRT(2.0)
      ROOT3=SQRT(3.0)
      AXIS=RCP(1)
      ALFA=RCP(4)*DEGRAD
      PSI=ARCOS(1.0/ROOT3)
      DIF=ARSIN(2.0*SIN(ALFA*0.5)/ROOT3)
      PHI=PSI-DIF
      PHIP=-COS(ALFA*0.5)/(ROOT3*COS(DIF))
      COSPHI=COS(PHI)
      SINPHI=SIN(PHI)
      BET=AXIS*COSPHI
      GAM=AXIS*SINPHI/ROOT2
      DEL=BET*(BET+GAM)-2.0*GAM**2
      ZDIA=(BET+GAM)/DEL
      ZOFF=-GAM/DEL
      SN(1)=COSPHI
      CS(1)=SINPHI/ROOT2
      SN(2)=-AXIS*SINPHI*PHIP
      CS(2)= AXIS*COSPHI*PHIP/ROOT2
      DO 1304 I = 1,6
      DO 1304 J = 1,6
 1304 DCELL(I,J)=0.0
      DO 1306 I = 1,2
      DELP=(2.0*BET+GAM)*SN(I)+(BET-4.0*GAM)*CS(I)
      ZDIAP=((SN(I)+CS(I))*DEL-(BET+GAM)*DELP)/DEL**2
      ZOFFP=(GAM*DELP-DEL*CS(I))/DEL**2
      DIA=0.5*ZDIA*ZDIAP+ZOFF*ZOFFP
      OFF=ZDIA*ZOFFP+ZOFF*(ZDIAP+ZOFFP)
      K=3*I-2
      DO 1305 J = 1,3
      DCELL(K,J  )=DIA
 1305 DCELL(K,J+3)=OFF
 1306 CONTINUE
      ASTAR=SQRT((BET+GAM)**2+2.0*GAM**2)/DEL
      ALFST=2.0*ATAN((ROOT3*COS(DIF))/1.0)/DEGRAD
      CELL=(BET-GAM)*DEL
      DO 1307 I = 1,3
      RECIP(I  )=ASTAR
      RECIP(I+3)=ALFST
      DO 1307 J = 1,3
      AT(I,J)=GAM
      AT(I,I)=BET
      ZT(I,J)=ZOFF
 1307 ZT(I,I)=ZDIA
      WRITE(IOO,52)
 1308 CONTINUE
      WRITE(IOO,51) (RCP(I),I=1,6),CELL,(RECIP(I),I=1,6),((AT(I,J),J=1,3
     1 ),(ZT(I,J),J=1,3),I=1,3)
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE MATINV (N,NFAIL)
      COMMON IOI,IOO,AM(861),V(41)
C     ********** SEGMENT 1 OF CHOLESKI INVERSION **********
C     ***** FACTOR MATRIX INTO LOWER TRIANGLE X TRANSPOSE ****
      K=1
      NFAIL=1
      IF(N.LT.1)  GO TO 300
      IF(N.GT.1)  GO TO 9
      AM(1)=1.0/AM(1)
      GO TO 204
C     ***** LOOP M OF A(L,M) *****
    9 DO 7 M=1,N
      IMAX=M-1
C     ***** LOOP L OF A(L,M) *****
      DO 6 L=M,N
      SUMA=0.0
      KLI=L
      KMI=M
      IF(IMAX.EQ.0)  GO TO 23
C     *****SUM OVER I=1,M-1 A(L,I)*A(M,I) *****
      DO 2 I=1,IMAX
      SUMA=SUMA+AM(KLI)*AM(KMI)
      J=N-I
      KLI=KLI+J
    2 KMI=KMI+J
C     *****TERM=C(L,M)-SUM *****
   23 TERM=AM(K)-SUMA
      IF(L.GT.M)  GO TO 5
      IF(TERM.GT.0.0)  GO TO 4
      NFAIL=K
      GO TO 300
C     ***** A(M,M)=SQRT(TERM) *****
    4 DENOM=SQRT (TERM)
      AM(K)=DENOM
      GO TO 6
C     ***** A(L,M)=TERM/A(M,M) *****
    5 AM(K)=TERM/DENOM
    6 K=K+1
    7 CONTINUE
C     ********** SEGMENT 2 OF CHOLESKI INVERSION **********
C     *****INVERSION OF TRIANGULAR MATRIX*****
  100 AM(1)=1.0/AM(1)
      KDM=1
C     ***** STEP L OF B(L,M) *****
      DO 104 L=2,N
      KDM=KDM+N-L+2
C     ***** RECIPROCAL OF DIAGONAL TERM *****
      TERM = 1.0/AM(KDM)
      AM(KDM)=TERM
      KMI=0
      KLI=L
      IMAX=L-1
C     ***** STEP M OF B(L,M) *****
      DO 103 M=1,IMAX
      K=KLI
C     ***** SUM TERMS *****
      SUMA=0.0
      DO 102 I=M,IMAX
      II=KMI+I
      SUMA=SUMA-AM(KLI)*AM(II)
  102 KLI=KLI+N-I
C     ***** MULT SUM * RECIP OF DIAGONAL *****
      AM(K)=SUMA*TERM
      J=N-M
      KLI=K+J
  103 KMI=KMI+J
  104 CONTINUE
C     ********** SEGMENT 3 OF CHOLESKI INVERSION **********
C     *****PREMULTIPLY LOWER TRIANGLE BY TRANSPOSE*****
  200 K=1
      DO 203 M=1,N
      KLI=K
      DO 202 L=M,N
      KMI=K
      IMAX=N-L+1
      SUMA=0.0
      DO 201 I=1,IMAX
      SUMA=SUMA+AM(KLI)*AM(KMI)
      KLI=KLI+1
  201 KMI=KMI+1
      AM(K)=SUMA
  202 K=K+1
  203 CONTINUE
  204 NFAIL=0
  300 RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE FINISH(IWHY)
      COMMON IOI,IOO
      IF(IWHY.NE.12) GO TO 3
      WRITE(IOO,2)
      STOP
    3 WRITE(IOO,1) IWHY
      STOP
    1 FORMAT(' STOPPED FOR REASON',I6/24X'1 : NO PARAMETERS TO BE VARIED
     1'/24X'2 : CONTRIBUTION OUTSIDE BLOCK ARRAY, REORDER REFLECTIONS'/2
     24X'3 : WRONG NUMBER OF PARAMETERS'/24X'4 : TOO MANY PARAMETERS'/24
     3X'5 : TOO MANY ATOMS'/24X'6 : TOO MANY REFLECTIONS'/24X'7 : TOO WI
     4DE A SCAN'/24X'8 : TOO MANY VARIABLES'/24X'9 : UNPHYSICAL TEMPERAT
     5URE FACTOR'/23X'10 : SINGULAR LEAST SQUARES MATRIX'/23X'11 : ZERO
     6DIAGONAL MATRIX ELEMENT'/23X'12 : SUCCESSFUL RUN'/23X'13 : MINPT =
     7 0'/)
    2 FORMAT(15H SUCCESSFUL RUN)
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE PLOTLP
      COMMON IOI,IOO,AM(861),V(41),DEGRAD,RCP(6),DCELL(6,6),AT(3,3),
     1 ZT(3,3),IRHOMB,LATTIS,MINRH,MAXRH,NBRAGG,INDEX(5,1500),STEP,
     2 WISE(3100,2),MINPT,MAXPT,PTMIN,PTMAX,RMAX,ISCREW(3),IGLIDE(3),
     3 NREFS,IC,NC,NV,NVM,NP,NM,ITF,ISTOP,KEY,NPC,NA,LVAR(41),TWOPIS,
     4 BETA(6,40),P(400),KI(400),XYZ(3,40),SF(40),AI(40),TO,PWISE(9),
     5 BH,DADAI(40),DADX(3,40),DADB(6,40),H(3),NF,ICENT,NS,NREFL,RHO,
     6 NOUT,NPM,TWOPI,YC,SINTHL,DC(400),FX(32,10),TS(3,24),IS(2,3,24),
     7 NO,NATOM,NNA,NNB,NNC,NND,NNE,NNF,NNG,NNH,LUM,NVL,YSUM,YMASK,
     8 YAV,NFRAC,NMAX,NEXT,IOF,IW,INHKL
      DIMENSION LABEL(12),IOUT(121),IDOTS(121)
      DATA ISPACE/1H /,IDOT/1H./,IPLUS/1HX/,ISTAR/1H*/,IOBS/1HO/,IM/1H-/
  100 FORMAT(3X,12I10)
  102 FORMAT(1X,121A1,I7)
  103 FORMAT(//' MULTIPLY (ABSCISSAE+PARAMETER 5) BY',F8.4,':(STEP)'//)
  104 FORMAT(1X,121A1,3I3)
  105 FORMAT(' REORDERED REFLECTIONS'/10(I7,2I3))
C     REORDER REFLECTIONS
      DO 399 NREFL = 2,NBRAGG
      INDX=INDEX(5,NREFL)
      K=NREFL
      J=0
      DO 394 I = 1,20
      IF(K.EQ.1) GO TO 395
      K=K-1
      IF(INDX.GE.INDEX(5,K)) GO TO 395
  394 J=J+1
  395 IF(J.EQ.0) GO TO 399
      DO 396 I = 1,5
  396 IOUT(I)=INDEX(I,NREFL)
      DO 397 I = 1,J
      K=NREFL-I
      L=K+1
      DO 397 M = 1,5
  397 INDEX(M,L)=INDEX(M,K)
      DO 398 I = 1,5
  398 INDEX(I,K)=IOUT(I)
  399 CONTINUE
      WRITE(6,105) ((INDEX(I,J),I=1,3),J=1,NBRAGG)
      WMAX=0.0
      DO 401 I = MINPT,MAXPT
      IF(WISE(I,1).GT.WMAX) WMAX=WISE(I,1)
      IF(WISE(I,2).GT.WMAX) WMAX=WISE(I,2)
  401 CONTINUE
      FSCALE=WMAX/116.
      DO 402 I = 1,12
  402 LABEL(I)=10.0*FSCALE*FLOAT(I)
      WRITE(6,103) STEP
      WRITE(6,100) LABEL
      DO 403 J = 1,121
      IDOTS(J)=ISPACE
      IF(MOD(J,10).EQ.1) IDOTS(J)=IDOT
  403 CONTINUE
      WRITE(6,102) IDOTS
      INDEX(5,NBRAGG+1)=NMAX+1
      DO 404 IFIND = 1,NBRAGG
      IF(INDEX(5,IFIND).GE.MINPT) GO TO 405
  404 CONTINUE
  405 DO 415 I = MINPT,MAXPT
      DO 406 J = 2,120
  406 IOUT(J)=ISPACE
      IOUT(1)=IDOT
      IOUT(121)=IDOT
      IF(MOD(I,50).NE.0) GO TO 408
      DO 407 J = 1,121
  407 IOUT(J)=IDOTS(J)
  408 MOBS=WISE(I,1)/FSCALE+1.5
      MCAL=WISE(I,2)/FSCALE+1.5
      IREFL=0
  409 IF(INDEX(5,IFIND).GT.I) GO TO 410
      IREFL=IREFL+1
      IFIND=IFIND+1
      GO TO 409
  410 IF(IREFL.EQ.0) GO TO 412
      IREFL=121-IREFL
      DO 411 J = IREFL,120
  411 IOUT(J)=IM
  412 IF(MOBS.GT.0) IOUT(MOBS)=IOBS
      IF(MCAL.GT.0) IOUT(MCAL)=ISTAR
      IF(MOBS.LT.0) IOUT(2)=IM
      IF(IOUT(MOBS).EQ.ISTAR) IOUT(MOBS)=IPLUS
      IF(IREFL.EQ.0) GO TO 413
C     PRINTS LAST INDEX SET FOR THIS POINT
      IFINDM=IFIND-1
      WRITE(6,104) IOUT,(INDEX(J,IFINDM),J=1,3)
      GO TO 415
  413 IF(MOD(I,10).EQ.0) GO TO 414
      WRITE(6,102) IOUT
      GO TO 415
  414 WRITE(6,102) IOUT,I
  415 CONTINUE
      WRITE(6,102) IDOTS
      WRITE(6,100) LABEL
      WRITE(6,103) STEP
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE PROD (A,B,C,N)
      DIMENSION A(3,40),B(3,3),C(3,40)
      DO 1 I=1,3
      DO 1 J=1,N
    1 A(I,J)=B(I,1)*C(1,J)+B(I,2)*C(2,J)+B(I,3)*C(3,J)
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE BORDER
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      FUNCTION ARSIN(X)
      ARSIN=ASIN(X)
      RETURN
      END
C
C----------------------------------------------------------------------
C
      FUNCTION ARCOS(X)
      ARCOS=ACOS(X)
      RETURN
      END
C
C----------------------------------------------------------------------
C
 
 
 
 
 
 
 
 
 
 
