c
c
c     FOR008 is the output unit for the calculated helix
c     FOR009 is input; including control cards and Diamond list
c     FOR010 is the program output
c
c
C0    HELIX-GENERATING PROGRAM FROM JMR, REVISED BY HRD(AUGUST 25,1980)
C0
C1    CARD 1: TITLE(20A4). (LEAVE LEFTMOST CHARACTER OF TITLE BLANK.)
C1
C2    CARD 2: FORMAT TO READ FRACTIONAL COORDINATES X,Y,Z AND ATOM NAME(20A4);
C2            IT IS EASIEST TO USE DIAMOND LIST FORMAT WITH CARTESIAN
C2            COORDINATES AND NAME(I)=IDW(3F10.5,11X,A4).
C2
C3    CARD 3: NO. OF INPUT ATOMS TO BE READ, AND FLAG FOR PRINTING INPUT
C3            COORDINATES(NO=0,YES=1): READ NAT,IFLAGP(2I5).
C3
C4    CARDS 4:ATOMIC COORDINATES X,Y,Z, AND ATOM NAME IN CARD 2 FORMAT.
C4
C5    CARD 5: CELL CONSTANTS A,B,C,ALPHA,BETA,GAMMA(6F10.5).
C5
C6    CARDS 6:ATOM PAIRS TO BE USED IN DETERMINING THE HELIX AXIS:
C6            READ N1,N2,NTERM,FIELD(3I5,15A4), WHERE N1,N2 ARE THE ATOM NOS.
C6            IN ORDER OF CARD 4 INPUT, NTERM=0 OR BLANK EXCEPT FOR THE
C6            LAST ATOM PAIR, AND FIELD(15) IS A COMMENT SPACE.
C6
C6    THE VECTORS ARE READ (N1,N2) WHERE N1,N2=ATOM NOS.:  THE VECTOR
C6    CALCULATED IS (N2-N1). BE SURE TO BE CONSISTENT IN YOUR CHOICE OF
C6    VECTOR DIRECTION.
C6
C7    CARD 7: NPMIN,NPMAX,IFPUN(3I5); NPMIN IS MINIMUM POWER OF THE HELIX
C7            OPERATOR, NPMAX IS MAXIMUM POWER OF THE HELIX OPERATOR, AND
C7            IFPUN.NE.0 FOR OUTPUT ON UNIT 'IPU'.
C7
C7    AS PRESENTLY WRITTEN IN THE 'PUTATM' SUBROUTINE, THE PROGRAM OUTPUTS
C7    ON UNIT 'IPU' FOR (IFPUN.NE.0): (IN DIAMOND LIST FORMAT)(FB,132,3300)
C7       1) THE TITLE RECORD('1'/1X,20A4/'0'),
C7       2) CARTESIAN COORDINATES X,Y,Z,NAME(I)=IDW,ENTRY
C7          (3F10.5,11X,A4,5X,I5,77X)
C7    FOR ALL EVEN POWERS OF THE HELIX OPERATOR AND ALL INPUT ATOMS.
C7
C    *****PRESENT ARRAY LIMITS ALLOW 496 INPUT ATOMS AND 50 HELIX VECTORS*****
C
      COMMON NAT,NVEC,ATM(3,500),NAME(500),X(3,50),VEC(3,50),
     *TITL(20),IN,IPR,IPU,PIFAC
      INTEGER LS1(50),LS2(50)
      REAL ROT(3,3),A(50,2),G(50),XYC(2),SGXYC(2),XYZC(3)
      REAL ROTOT(3,3),WORK(3,3),OMAT(3,3),OUHV(3)
      REAL X1(3),X2(3),UHV(3),XLS(50,50),ELS(50),UXYZ(3,3)
      DATA UXYZ/1.0, 3*0.0, 1.0, 3*0.0, 1.0/
      PIFAC=3.1415926/180.0
C
C     I/O UNIT ASSIGNMENTS:CARD READ,PRINTER,AND IFPUN.NE.0
C
      IN=9
      IPR=10
      IPU=8
C
      CALL GETATM(OMAT)
      CALL GETHVC
      CALL MLSP(NVEC,VEC,UHV,DZ,SIGPL)
      TOL=1.0E-5
      IF((ABS(UHV(1)).LE.TOL).AND.(ABS(UHV(2)).LE.TOL)) GO TO 200
      ANG=ACOS(UHV(1)/SQRT(UHV(1)**2+UHV(2)**2))
      IF(UHV(2).LE.0.0) ANG=-ANG
      CALL DEFROT(UXYZ(1,3),ANG,ROT)
      CALL MATMUL(ROTOT,ROT,UXYZ,WORK,3,3,3)
      DO 100 I=1,NAT
      CALL ROTATE(ROT,ATM(1,I))
  100 CONTINUE
      DO 110 I=1,NVEC
      CALL ROTATE(ROT,X(1,I))
      CALL ROTATE(ROT,VEC(1,I))
  110 CONTINUE
  200 CONTINUE
      ANG=ACOS(UHV(3))
      CALL DEFROT(UXYZ(1,2),ANG,ROT)
      CALL MATMUL(ROTOT,ROT,ROTOT,WORK,3,3,3)
      DO 120 I=1,NAT
      CALL ROTATE(ROT,ATM(1,I))
  120 CONTINUE
      DO 130 I=1,NVEC
      CALL ROTATE(ROT,X(1,I))
      CALL ROTATE(ROT,VEC(1,I))
  130 CONTINUE
      DO 210 I=1,NVEC
      A(I,1)=VEC(1,I)*2.0
      A(I,2)=VEC(2,I)*2.0
      G(I)=(VEC(1,I) + X(1,I))**2 + (VEC(2,I) + X(2,I))**2
      G(I)=G(I)-X(1,I)**2 -X(2,I)**2
  210 CONTINUE
      CALL LLSQ(NVEC,2,50,2,G,A,XYC,SGXYC,NSING,XLS,ELS,LS1,LS2)
      IF(NSING.EQ.0) GO TO 220
      WRITE(IPR,2200)
 2200 FORMAT(' >>>>>>>>>>> CENTER MATRIX SINGULAR')
      STOP 16
  220 CONTINUE
      XYZC(1)=XYC(1)
      XYZC(2)=XYC(2)
      XYZC(3)=0.0
      CALL MATMUL(XYZC,XYZC,ROTOT,WORK,1,3,3)
      WRITE(IPR,2220)
 2220 FORMAT(' IN ORTHONORMAL COORDINATES, ')
      WRITE(IPR,2210) (UHV(JI),XYZC(JI),JI =1,3)
 2210 FORMAT(' HELIX AXIS IN PARAMETRIC FORM:'/' X = ',F12.5,'*S + ',
     *F12.5/' Y = ',F12.5,'*S + ',F12.5/' Z = ',F12.5, '*S + ',F12.5/)
      CALL MATMUL(XYZC,OMAT,XYZC,WORK,3,3,1)
      CALL MATMUL(OUHV,OMAT,UHV,WORK,3,3,1)
      WRITE(IPR,2230)
 2230 FORMAT(' IN ORIGINAL CRYSTAL COORDINATES,')
      WRITE(IPR,2210) (OUHV(JI),XYZC(JI),JI=1,3)
      DO 250 I=1,2
      DO 230 J=1,NAT
      ATM(I,J)=ATM(I,J)-XYC(I)
  230 CONTINUE
      DO 240 J=1,NVEC
      X(I,J)=X(I,J)-XYC(I)
  240 CONTINUE
  250 CONTINUE
      THETA=0.0
      SGTHET=0.0
      DO 310 I=1,NVEC
      DO 300 J=1,3
      X1(J)=X(J,I)
      X2(J)=X(J,I)+VEC(J,I)
  300 CONTINUE
      CALL POLAR(X1)
      CALL POLAR(X2)
      T21=X2(2)-X1(2)
      T21P=T21+PIFAC*360.0
      T21M=T21-PIFAC*360.0
      IF(ABS(T21P).LT.ABS(T21)) T21=T21P
      IF(ABS(T21M).LT.ABS(T21)) T21=T21M
      THETA=THETA+T21
      SGTHET=SGTHET+T21*T21
  310 CONTINUE
      THETA=THETA/FLOAT(NVEC)
      SGTHET=(SGTHET-FLOAT(NVEC)*THETA*THETA)/FLOAT(NVEC-1)
      SGTHET=SQRT(SGTHET)/PIFAC
      SIG=0.0
      IF(NVEC.LE.3) GO TO 340
      DO 330 I=1,NVEC
      DO 320 J=1,3
      X1(J)=X(J,I)
      X2(J)=X(J,I)+VEC(J,I)
  320 CONTINUE
      CALL POLAR(X1)
      CALL POLAR(X2)
      SIG=SIG+(X1(3)+DZ-X2(3))**2 + X1(1)**2 + X2(1)**2
     *-2.0*X1(1)*X2(1)*COS(X1(2)+THETA-X2(2))
  330 CONTINUE
      SIG=SQRT(SIG/FLOAT(NVEC-3))
  340 CONTINUE
      IF(DZ.LT.0.0) THETA=-THETA
      IF(DZ.LT.0.0) DZ=-DZ
      ANG=THETA/PIFAC
      WRITE(IPR,2300) ANG,DZ,SIG
 2300 FORMAT(' >>>>>> HELIX ROTATION: ',F8.3,5X,'DISPLACEMENT: ',
     *F8.4/' ',10X,'STATISTICS:'/' ',20X,'OVERALL STANDARD DEV.:  ',
     *F8.4)
      WRITE(IPR,2310) SGXYC,SIGPL,SGTHET
 2310 FORMAT(21X,'SIGMA(X): ',F7.4,', SIGMA(Y): ',F7.4,
     *', SIGMA(Z)=SIGMA(DISPLACEMENT): ',F7.4,', SIGMA(ROTATION): ',
     *F7.3)
      PTURN=ABS(360.0/ANG)
      WRITE(IPR,2320)PTURN
 2320 FORMAT(11X,'THERE ARE ',F5.2,' RESIDUES PER TURN'/)
      CALL PUTATM(THETA,DZ)
      WRITE(IPR,2400)
 2400 FORMAT(' NORMAL END OF JOB: YOU HAVE GIVEN BIRTH TO A HELIX')
      STOP
      END
      SUBROUTINE GETATM(OMAT)
C
C     READ IN PARAMETERS AND ATOMIC COORDINATES
C
      COMMON NAT,NVEC,ATM(3,500),NAME(500),X(3,50),VEC(3,50),
     *TITL(20),IN,IPR,IPU,PIFAC
      DIMENSION FMT(20),W1(3),W2(3),CELL(6),OMAT(3,3),OATM(3,500)
      DATA NORIG/'ORGN'/,NOX/'A AX'/,NOY/'B AX'/,NOZ/'C AX'/
      READ(IN,1000)TITL
 1000 FORMAT(20A4)
      WRITE(IPR,2000) TITL
 2000 FORMAT(20A4/)
      WRITE(IPR,2222)
 2222 FORMAT(' INPUT FRAC.COORDS. X,Y,Z AND ATOM NAME IN FORMAT:')
      READ(IN,1000)FMT
      WRITE(IPR,2000)FMT
      READ(IN,1010) NAT,IFLAGP
 1010 FORMAT(2I5)
      WRITE(IPR,2010) NAT,IFLAGP
 2010 FORMAT(' WILL READ ',I5,' ATOMS: PRINT FLAG=',I5/)
C
C     >>>>>>>>SET MAX. NO. OF ATOMS=ARRAY SIZE-4<<<<<<<<<<<<<<<<<
      IF(NAT.LE.496) GO TO 100
C
      WRITE(IPR,2100)
 2100 FORMAT(' THAT EXCEEDS THE ATM(I,J),NAME(I),AND OATM(I,J) ARRAYS')
      STOP 16
  100 CONTINUE
      DO 10 I=1,NAT
C
C  THE FOLLOWING READ STATEMENT READS IN COORDINATES X,Y,Z AND IDENTI-
C   CATION NUMBER.  IF FOR ANY REASON YOU MUST PERMUTE COORDINATES,
C   THIS CAN BE DONE JUST AFTER THE READ STATEMENT.
C   IF YOUR INPUT FILE HAS THE IDENTIFICATION NUMBER BEFORE THE
C   X,Y,Z COORDINATES, THEN REWRITE THE FOLLOWING READ STATEMENT,
C   WITH THE FORMAT LISTING AT THE BEGINNING OF INPUT FILE 9 ADJUSTED
C   ACCORDINGLY.
C
      READ(IN,FMT) (ATM(J,I),J=1,3),NAME(I)
      IF(IFLAGP.NE.1) GO TO 10
      WRITE(IPR,2015) I,(ATM(J,I),J=1,3),NAME(I)
 2015 FORMAT(1X,I5,5X,3(F9.4,1X),1X,A4)
   10 CONTINUE
      ATM(1,NAT+1)=0.0
      ATM(2,NAT+1)=0.0
      ATM(3,NAT+1)=0.0
      ATM(1,NAT+2)=1.0
      ATM(2,NAT+2)=0.0
      ATM(3,NAT+2)=0.0
      ATM(1,NAT+3)=0.0
      ATM(2,NAT+3)=1.0
      ATM(3,NAT+3)=0.0
      ATM(1,NAT+4)=0.0
      ATM(2,NAT+4)=0.0
      ATM(3,NAT+4)=1.0
      NAME(NAT+1)=NORIG
      NAME(NAT+2)=NOX
      NAME(NAT+3)=NOY
      NAME(NAT+4)=NOZ
      NAT=NAT+4
      READ(IN,1020) CELL
 1020 FORMAT(6F10.5)
      WRITE(IPR,2333)
 2333 FORMAT(/' CELL CONSTANTS A,B,C,ALPHA,BETA,GAMMA ARE:')
      WRITE(IPR,2020) CELL
 2020 FORMAT(1X,6(F7.3,1X)/)
      CALL GETMAT(OMAT,CELL)
      WRITE(IPR,2025) IFLAGP
 2025 FORMAT(' COORDINATES IN ORTHONORMAL SYSTEM: PRINT FLAG=',I5/)
      DO 20 I=1,NAT
      CALL ROTATE(OMAT,ATM(1,I))
      IF(IFLAGP.NE.1) GO TO 36
      WRITE(IPR,2015) I,(ATM(J,I),J=1,3),NAME(I)
   36 CONTINUE
      DO 20 J=1,3
      OATM(J,I)=ATM(J,I)
   20 CONTINUE
      CALL MINV(OMAT,3,DETOMT,W1,W2)
      RETURN
      END
      SUBROUTINE GETHVC
C
C     GET VECTORS RELATING EQUIVALENT ATOMS
C
      COMMON NAT,NVEC,ATM(3,500),NAME(500),X(3,50),VEC(3,50),
     *TITL(20),IN,IPR,IPU,PIFAC
      DIMENSION FIELD(15)
      WRITE(IPR,2444)
 2444 FORMAT(' **************************************'/)
      WRITE(IPR,2000)
 2000 FORMAT(' ATOM PAIRS USED TO DETERMINE THE HELIX AXIS ARE:')
      NVEC=1
   10 READ(IN,1000) N1,N2,NTERM,FIELD
 1000 FORMAT(3I5,15A4)
      WRITE(IPR,2010) N1,NAME(N1),N2,NAME(N2),FIELD
 2010 FORMAT(1X,2(I3,1X,A4,8X),15A4)
      DO 20 I=1,3
      X(I,NVEC)=ATM(I,N1)
      VEC(I,NVEC)=ATM(I,N2)-ATM(I,N1)
   20 CONTINUE
      IF(NTERM.NE.0) GO TO 100
      NVEC=NVEC+1
      IF(NVEC.LE.50) GO TO 10
      WRITE(IPR,2020)
 2020 FORMAT(' ARRAYS EXCEEDED BY INPUT OF GT.50 VECTORS')
      STOP 16
  100 WRITE(IPR,2444)
      IF(NVEC.GE.3) RETURN
      WRITE(IPR,2030)
 2030 FORMAT(' LESS THAN 3 VECTORS: HELIX AXIS CANNOT BE DETERMINED')
      STOP 16
      END
      SUBROUTINE PUTATM(ANG,DZ)
C
C     OUTPUT ATOMIC POSITIONS FOR VARIOUS POWERS OF THE HELIX OPERATOR
C
      COMMON NAT,NVEC,ATM(3,500),NAME(500),X(3,50),VEC(3,50),
     *TITL(20),IN,IPR,IPU,PIFAC
      DIMENSION XA(3)
      NPMIN=0
      NPMAX=0
      IFPUN=0
      DO 5 I=1,NAT
      CALL POLAR(ATM(1,I))
    5 CONTINUE
      READ(IN,1000)NPMIN,NPMAX,IFPUN
 1000 FORMAT(3I5)
      NPMAX=NPMAX-NPMIN+1
C
C     OUTPUT ON 'IPR' ALL POWERS OF THE HELIX OPERATOR, BOTH ODD AND EVEN
C
      DO 10 I=1,NPMAX
      POWER=FLOAT(I+NPMIN-1)
      IP=I+NPMIN-1
      AD=ANG*POWER
      ZD=DZ*POWER
      WRITE(IPR,2444)
 2444 FORMAT(' *****************************************************'/)
      WRITE(IPR,2000) TITL,IP
 2000 FORMAT(1X,20A4/'   NO.',5X,'    R',6X,'THETA',7X,'Z',5X,'NAME  ',
     *4X,'X',8X,'Y',8X,'Z',8X,'POWER = ',I3/)
      DO 10 J=1,NAT
      R=ATM(1,J)
      T=ATM(2,J)+AD
      Z=ATM(3,J)+ZD
      XC=R*COS(T)
      Y=R*SIN(T)
      T=PRIANG(T/PIFAC)
      WRITE(IPR,2010) J,R,T,Z,NAME(J),XC,Y,Z
 2010 FORMAT(1X,I5,5X,3(F8.3,1X),1X,A4,2X,3(F8.3,1X))
   10 CONTINUE
C
C     OUTPUT ON UNIT 'IPU':ONLY EVEN POWERS OF THE HELIX OPERATOR IN
C    *DIAMOND LIST FORMAT(RECORD LENGTH=132)
C
      IF(IFPUN.EQ.0) GO TO 70
      NOUT=0
      WRITE(IPR,2666)
 2666 FORMAT(' THE FOLLOWING DIAMOND LIST HAS BEEN OUTPUT ON UNIT IPU:')
      WRITE(IPR,2555)TITL
      WRITE(IPU,2555)TITL
 2555 FORMAT('0'/20A4/'0')
      DO 20 I=1,NPMAX
      POWER=FLOAT(I+NPMIN-1)
      IP=I+NPMIN-1
      IPMOD=MOD(IP,2)
      IF(IPMOD.NE.0) GO TO 20
      AD=ANG*POWER
      ZD=DZ*POWER
      NAT2=NAT-4
      DO 20 J=1,NAT2
      R=ATM(1,J)
      T=ATM(2,J)+AD
      Z=ATM(3,J)+ZD
      XC=R*COS(T)
      Y=R*SIN(T)
      T=PRIANG(T/PIFAC)
      NOUT=NOUT+1
      WRITE(IPR,2020) XC,Y,Z,NAME(J),NOUT
      WRITE(IPU,2020) XC,Y,Z,NAME(J),NOUT
 2020 FORMAT(3F10.5,11X,A4,5X,I5,77X)
   20 CONTINUE
   70 CONTINUE
      RETURN
      END
      FUNCTION PRIANG(ANGLE)
C
C     SET ANGLE T TO BETWEEN 0.0 AND 360.0
C
      T=ANGLE
   30 IF(T.LT.360.0) GO TO 40
      T=T-360.0
      GO TO 30
   40 IF(T.GE.0.0) GO TO 50
      T=T+360.0
      GO TO 40
   50 CONTINUE
      PRIANG=T
      RETURN
      END
      SUBROUTINE MLSP(NPOINT,X,U,D,SIG)
C
C     GENERAL LEAST-SQUARES PLANE ROUTINE
C
      DIMENSION X(3,NPOINT),U(3),XAV(3),XM(3,3),R(3,3),YM(3,3),ZM(9)
      EQUIVALENCE (XM(1,1),ZM(1))
      IF(NPOINT.LT.3) STOP 16
      DO 10 I=1,3
      XAV(I)=0.0
   10 CONTINUE
      DO 20 I=1,3
      DO 20 J=1,NPOINT
      XAV(I)=XAV(I)+X(I,J)
   20 CONTINUE
      DO 30 I=1,3
      XAV(I)=XAV(I)/FLOAT(NPOINT)
   30 CONTINUE
      DO 40 I=1,3
      DO 40 J=1,3
      XM(I,J)=0.0
      YM(I,J)=0.0
   40 CONTINUE
      DO 50 I=1,3
      DO 50 J=1,3
      DO 50 K=1,NPOINT
      YM(I,J)=YM(I,J)+(X(I,K)-XAV(I))*(X(J,K)-XAV(J))
   50 CONTINUE
      ZM(1)=YM(1,1)
      ZM(2)=YM(1,2)
      ZM(3)=YM(2,2)
      ZM(4)=YM(1,3)
      ZM(5)=YM(2,3)
      ZM(6)=YM(3,3)
      CALL EIGEN(XM,R,3,0)
      EVAL=ZM(6)
      SIG=0.0
      IF(NPOINT.LE.3) GO TO 60
      IF(EVAL.LT.0.0) GO TO 60
      SIG=SQRT(EVAL/FLOAT(NPOINT-3))
   60 CONTINUE
      DO 70 I=1,3
      U(I)=R(I,3)
   70 CONTINUE
      UMAG=DOTP(U,U)
      DO 80 I=1,3
      U(I)=U(I)/UMAG
   80 CONTINUE
      D=0.0
      DO 90 I=1,3
      D=D+U(I)*XAV(I)
   90 CONTINUE
      RETURN
      END
      SUBROUTINE DEFROT(VECTOR,ANGLE,ROTOR)
C
C     DEFINE A GENERAL ROTATION MATRIX
C
      DIMENSION VECTOR(3),DIRCOS(3),ROTOR(3,3)
      DO 10 I=1,3
      DO 10 J=1,3
   10 ROTOR(I,J)=0.0
      C=COS(ANGLE)
      S=SIN(ANGLE)
      VMAG=SQRT(DOTP(VECTOR,VECTOR))
      DO 20 I=1,3
   20 DIRCOS(I)=VECTOR(I)/VMAG
      ROTOR(1,2)=DIRCOS(3)*S
      ROTOR(2,1)=-ROTOR(1,2)
      ROTOR(3,1)=DIRCOS(2)*S
      ROTOR(1,3)=-ROTOR(3,1)
      ROTOR(2,3)=DIRCOS(1)*S
      ROTOR(3,2)=-ROTOR(2,3)
      DO 30 I=1,3
      DO 30 J=1,3
   30 ROTOR(I,J)=ROTOR(I,J)+DIRCOS(I)*DIRCOS(J)*(1.0-C)
      DO 40 I=1,3
   40 ROTOR(I,I)=ROTOR(I,I)+C
      RETURN
      END
      SUBROUTINE ROTATE(ROTOR,VECTOR)
C
C     ROTATE A VECTOR VIA THE MATRIX ROTOR
C
      DIMENSION ROTOR(3,3),VECTOR(3),WORK(3)
      DO 10 I=1,3
      WORK(I)=0.0
   10 CONTINUE
      DO 20 I=1,3
      DO 20 J=1,3
      WORK(I)=WORK(I)+ROTOR(I,J)*VECTOR(J)
   20 CONTINUE
      DO 30 I=1,3
      VECTOR(I)=WORK(I)
   30 CONTINUE
      RETURN
      END
      SUBROUTINE LLSQ(NOBS,NPARM,NOD,NPD,G,A,P,SIGP,NSING,M,E,L1,L2)
C
C     GENERAL LINEAR LEAST-SQUARES ROUTINE
C
      DIMENSION G(NOD),A(NOD,NPD),P(NPD),SIGP(NPD)
      REAL   M(NPARM,NPARM),E(NPARM),L1(NPARM),L2(NPARM)
      IF((NOBS.GT.NOD).OR.(NPARM.GT.NPD)) STOP 16
      DO 10 I=1,NPARM
      E(I)=0.0
      DO 10 J=1,NPARM
      M(I,J)=0.0
   10 CONTINUE
      DO 20 I=1,NPARM
      DO 20 J=1,NOBS
      E(I)=E(I)+A(J,I)*G(J)
      DO 20 K=1,NPARM
      M(I,K)=M(I,K)+A(J,I)*A(J,K)
   20 CONTINUE
      CALL MINV(M,NPARM,DET,L1,L2)
      NSING=0
      IF(DET.NE.0.0) GO TO 30
      NSING=1
      RETURN
   30 CONTINUE
      DO 40 I=1,NPARM
      P(I)=0.0
   40 CONTINUE
      DO 50 I=1,NPARM
      DO 50 J=1,NPARM
      P(I)=P(I)+M(I,J)*E(J)
   50 CONTINUE
      SIGSQ=0.0
      DO 70 I=1,NOBS
      DEL=G(I)
      DO 60 J=1,NPARM
      DEL=DEL-A(I,J)*P(J)
   60 CONTINUE
      SIGSQ=SIGSQ+DEL*DEL
   70 CONTINUE
      DO 80 I=1,NPARM
      SIGP(I)=SQRT(M(I,I)*SIGSQ)
   80 CONTINUE
      RETURN
      END
      SUBROUTINE POLAR(V)
C
C     CONVERT A CARTESIAN VECTOR V INTO CYLINDRICAL POLAR
C
      DIMENSION V(3)
      X=V(1)
      Y=V(2)
      V(1)=SQRT(X*X+Y*Y)
      V(2)=ACOS(X/V(1))
      IF(Y.LT.0.0) V(2)=-V(2)
      RETURN
      END
      SUBROUTINE GETMAT(BMAT,CELCON)
C
C     GET THE MATRIX TO TRANSFORM INTO CARTESIAN COORDINATES
C
      REAL   BMAT(3,3),CELCON(6),ANG(3)
      ANG(1)=CELCON(4)
      ANG(2)=CELCON(5)
      ANG(3)=CELCON(6)
      PIFAZE=3.1415926/180.0
      DO 10 I=1,3
   10 ANG(I)=ANG(I)*PIFAZE
      A1=CELCON(1)
      A2=CELCON(2)
      A3=CELCON(3)
      TEMP1=(COS(ANG(2))-COS(ANG(1))*COS(ANG(3)))
      TEMP2=SIN(ANG(1))*SIN(ANG(3))
      OMG=ACOS(TEMP1/TEMP2)
      BMAT(1,1)=A1*SIN(ANG(3))*SIN(OMG)
      BMAT(1,2)=0.0
      BMAT(1,3)=0.0
      BMAT(2,1)=A1*COS(ANG(3))
      BMAT(2,2)=A2
      BMAT(2,3)=A3*COS(ANG(1))
      BMAT(3,1)=A1*SIN(ANG(3))*COS(OMG)
      BMAT(3,2)=0.0
      BMAT(3,3)=A3*SIN(ANG(1))
      RETURN
      END
      SUBROUTINE EIGEN(A,R,N,MV)
C
C     COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC MATRIX
C
      DIMENSION A(1),R(1)
      RANGE=1.0E-6
      IF(MV-1) 10,25,10
   10 IQ=-N
      DO 20 J=1,N
      IQ=IQ+N
      DO 20 I=1,N
      IJ=IQ+I
      R(IJ)=0.0
      IF(I-J) 20,15,20
   15 R(IJ)=1.0
   20 CONTINUE
   25 ANORM=0.0
      DO 35 I=1,N
      DO 35 J=I,N
      IF(I-J) 30,35,30
   30 IA=I+(J*J-J)/2
      ANORM=ANORM+A(IA)*A(IA)
   35 CONTINUE
      IF(ANORM) 165,165,40
   40 ANORM=1.414*SQRT(ANORM)
      ANRMX=ANORM*RANGE/FLOAT(N)
      IND=0
      THR=ANORM
   45 THR=THR/FLOAT(N)
   50 L=1
   55 M=L+1
   60 MQ=(M*M-M)/2
      LQ=(L*L-L)/2
      LM=L+MQ
   62 IF(ABS(A(LM))-THR) 130,65,65
   65 IND=1
      LL=L+LQ
      MM=M+MQ
      X=0.5*(A(LL)-A(MM))
   68 Y=-A(LM)/SQRT(A(LM)*A(LM)+X*X)
      IF(X) 70,75,75
   70 Y=-Y
   75 SINX=Y/SQRT(2.0*(1.0+(SQRT(1.0-Y*Y))))
      SINX2=SINX*SINX
   78 COSX=SQRT(1.0-SINX2)
      COSX2=COSX*COSX
      SINCS=SINX*COSX
      ILQ=N*(L-1)
      IMQ=N*(M-1)
      DO 125 I=1,N
      IQ=(I*I-I)/2
      IF(I-L) 80,115,80
   80 IF(I-M) 85,115,90
   85 IM=I+MQ
      GO TO 95
   90 IM=M+IQ
   95 IF(I-L) 100,105,105
  100 IL=I+LQ
      GO TO 110
  105 IL=L+IQ
  110 X=A(IL)*COSX-A(IM)*SINX
      A(IM)=A(IL)*SINX+A(IM)*COSX
      A(IL)=X
  115 IF(MV-1) 120,125,120
  120 ILR=ILQ+I
      IMR=IMQ+I
      X=R(ILR)*COSX-R(IMR)*SINX
      R(IMR)=R(ILR)*SINX+R(IMR)*COSX
      R(ILR)=X
  125 CONTINUE
      X=2.0*A(LM)*SINCS
      Y=A(LL)*COSX2+A(MM)*SINX2-X
      X=A(LL)*SINX2+A(MM)*COSX2+X
      A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2)
      A(LL)=Y
      A(MM)=X
  130 IF(M-N) 135,140,135
  135 M=M+1
      GO TO 60
  140 IF(L-(N-1)) 145,150,145
  145 L=L+1
      GO TO 55
  150 IF(IND-1) 160,155,160
  155 IND=0
      GO TO 50
  160 IF(THR-ANRMX) 165,165,45
  165 IQ=-N
      DO 185 I=1,N
      IQ=IQ+N
      LL=I+(I*I-I)/2
      JQ=N*(I-2)
      DO 185 J=I,N
      JQ=JQ+N
      MM=J+(J*J-J)/2
      IF(A(LL)-A(MM)) 170,185,185
  170 X=A(LL)
      A(LL)=A(MM)
      A(MM)=X
      IF(MV-1) 175,185,175
  175 DO 180 K=1,N
      ILR=IQ+K
      IMR=JQ+K
      X=R(ILR)
      R(ILR)=R(IMR)
  180 R(IMR)=X
  185 CONTINUE

      RETURN
      END
      FUNCTION DOTP(X,Y)
C
C     TAKE THE DOT PRODUCT OF TWO CARTESIAN VECTORS
C
      DIMENSION X(3),Y(3)
      DOTP=0.0
      DO 10 I=1,3
      DOTP=DOTP+X(I)*Y(I)
   10 CONTINUE
      RETURN
      END
      SUBROUTINE MINV(A,N,D,L,M)
C
C     INVERT A MATRIX
C
      DIMENSION A(1),L(1),M(1)
      D=1.0
      NK=-N
      DO 80 K=1,N
      NK=NK+N
      L(K)=K
      M(K)=K
      KK=NK+K
      BIGA=A(KK)
      DO 20 J=K,N
      IZ=N*(J-1)
      DO 20 I=K,N
      IJ=IZ+I
   10 IF(ABS(BIGA)-ABS(A(IJ))) 15,20,20
   15 BIGA=A(IJ)
      L(K)=I
      M(K)=J
   20 CONTINUE
      J=L(K)
      IF(J-K) 35,35,25
   25 KI=K-N
      DO 30 I=1,N
      KI=KI+N
      HOLD=-A(KI)
      JI=KI-K+J
      A(KI)=A(JI)
   30 A(JI)=HOLD
   35 I=M(K)
      IF(I-K) 45,45,38
   38 JP=N*(I-1)
      DO 40 J=1,N
      JK=NK+J
      JI=JP+J
      HOLD=-A(JK)
      A(JK)=A(JI)
   40 A(JI)=HOLD
   45 IF(BIGA) 48,46,48
   46 D=0.0
      RETURN
   48 DO 55 I=1,N
      IF(I-K) 50,55,50
   50 IK=NK+I
      A(IK)=A(IK)/(-BIGA)
   55 CONTINUE
      DO 65 I=1,N
      IK=NK+I
      HOLD=A(IK)
      IJ=I-N
      DO 65 J=1,N
      IJ=IJ+N
      IF(I-K) 60,65,60
   60 IF(J-K) 62,65,62
   62 KJ=IJ-I+K
      A(IJ)=HOLD*A(KJ)+A(IJ)
   65 CONTINUE
      KJ=K-N
      DO 75 J=1,N
      KJ=KJ+N
      IF(J-K) 70,75,70
   70 A(KJ)=A(KJ)/BIGA
   75 CONTINUE
      D=D*BIGA
      A(KK)=1.0/BIGA
   80 CONTINUE
      K=N
  100 K=(K-1)
      IF(K) 150,150,105
  105 I=L(K)
      IF(I-K) 120,120,108
  108 JQ=N*(K-1)
      JR=N*(I-1)
      DO 110 J=1,N
      JK=JQ+J
      HOLD=A(JK)
      JI=JR+J
      A(JK)=-A(JI)
  110 A(JI)=HOLD
  120 J=M(K)
      IF(J-K) 100,100,125
  125 KI=K-N
      DO 130 I=1,N
      KI=KI+N
      HOLD=A(KI)
      JI=KI-K+J
      A(KI)=-A(JI)
  130 A(JI)=HOLD
      GO TO 100
  150 RETURN
      END
      SUBROUTINE MATMUL(PROD,XM,YM,WM,N,M,L)
C
C     MATRIX MULTIPLICATION
C
      DIMENSION XM(N,M),YM(M,L),PROD(N,L),WM(N,L)
      DO 10 I=1,N
      DO 10 J=1,L
      WM(I,J)=0.0
   10 CONTINUE
      DO 20 I=1,N
      DO 20 K=1,L
      DO 20 J=1,M
      WM(I,K)=WM(I,K)+XM(I,J)*YM(J,K)
   20 CONTINUE
      DO 30 I=1,N
      DO 30 J=1,L
      PROD(I,J)=WM(I,J)
   30 CONTINUE
      RETURN
      END