      SUBROUTINE DIRECT(SM,V,IPH,N)
      include 'param.inc'
      COMMON/F1/SMM(MSZ,MSZ),VX(MSZ)
      COMMON/PARAMS/GLB(16),AGLB(16),LGLB(16),PAR(NUPH,20),
     *APAR(NUPH,20),LPAR(NUPH,20),XL(NATS,11),LP(NATS,11),A(NATS,11),
     *PTR(NATS),RELAX(4),RATIO(2)
      CHARACTER*4 ATEXT(NATS),NTYP(NATS)
      COMMON/PARAC/ATEXT,NTYP
      DIMENSION V(6),SM(6,6)
      DO 1 J=1,6
      V(J)=PAR(IPH,J+5)
      K=LPAR(IPH,J+5)
      X=APAR(IPH,J+5)
      IF(K.EQ.0)X=0.
      DO 1 L=J,6
      M=LPAR(IPH,L+5)
      SM(L,J)=0.
      IF((M.EQ.0).OR.(K.EQ.0))GOTO 1
      SM(L,J)=SMM(M,K)*X*APAR(IPH,L+5)
    1 SM(J,L)=SM(L,J)
      CALL ESD(SM,V,1.)
      RETURN
      END
      SUBROUTINE ESD(SM,V,SUM)
      include 'param.inc'
      DIMENSION SM(6,6),V(6),S(6),R(6),E(6)
      I=0
      DO 1 IA=1,3
      I=I+1
      IB=IA+1
      IF(IB.GT.3)IB=IB-3
      IC=IB+1
      IF(IC.GT.3)IC=IC-3
      ID=IA+3
      IE=IB+3
      IP=IC+3
      FNUM=4.*V(IB)*V(IC)-V(ID)*V(ID)
      DEN=4.*V(IA)*V(IB)*V(IC)-V(IA)*V(ID)*V(ID)-V(IB)*V(IE)*V(IE)-V(IC)
     1*V(IP)*V(IP)+V(ID)*V(IE)*V(IP)
      S(I)=SQRT(FNUM/DEN)
      R(IA)=-S(I)**4
      DENL=1./(DEN*DEN)
      R(IB)=DENL*(4.*V(IC)*DEN-(4.*V(IA)*V(IC)-V(IE)*V(IE))*FNUM)
      R(IC)=DENL*(4.*V(IB)*DEN-(4.*V(IA)*V(IB)-V(IP)*V(IP))*FNUM)
      R(ID)=-DENL*(2.*V(ID)*DEN-(2.*V(IA)*V(ID)-V(IE)*V(IP))*FNUM)
      R(IE)=DENL*FNUM*(2.*V(IB)*V(IE)-V(ID)*V(IP))
      R(IP)=DENL*FNUM*(2.*V(IC)*V(IP)-V(ID)*V(IE))
      E(I)=ERROR(SM(1,1),R(1),SUM)
      E(I)=E(I)/(2.*S(I))
      FNUM=V(IE)*V(IP)-2.*V(IA)*V(ID)
      DEN=16.*V(IA)*V(IA)*V(IB)*V(IC)+V(IE)*V(IE)*V(IP)*V(IP)-4.*(V(IA)*
     1V(IB)*V(IE)*V(IE)+V(IC)*V(IA)*V(IP)*V(IP))
      DENL=1./(DEN*DEN)
      FNUML=FNUM*FNUM
      S(I+3)=FNUML/DEN
      R(IA)=-4.*FNUM*DENL*(V(ID)*DEN+FNUM*(8.*V(IA)*V(IB)*V(IC)-V(IB)*V(
     1IE)*V(IE)-V(IC)*V(IP)*V(IP)))
      R(IB)=-4.*FNUML*DENL*(4.*V(IA)*V(IA)*V(IC)-V(IA)*V(IE)*V(IE))
      R(IC)=-4.*FNUML*DENL*(4.*V(IA)*V(IA)*V(IB)-V(IA)*V(IP)*V(IP))
      R(ID)=-4.*FNUM/DEN*V(IA)
      R(IE)=2.*FNUM*DENL*(V(IP)*DEN-FNUM*(V(IE)*V(IP)*V(IP)-4.*V(IA)*V(I
     1B)*V(IE)))
      R(IP)=2.*FNUM*DENL*(V(IE)*DEN-FNUM*(V(IE)*V(IE)*V(IP)-4.*V(IC)*V(I
     1A)*V(IP)))
      E(I+3)=ERROR(SM(1,1),R(1),SUM)
      IF(S(I+3).EQ.0.0) GOTO 3
      E(I+3)=E(I+3)/(2.*SQRT(S(I+3)*(1.-S(I+3))))*180./3.14159265359
      S(I+3)=ATAN(SQRT((1.-S(I+3))/S(I+3)))*180./3.14159265359
      GOTO 4
    3 S(I+3)=90.
      E(I+3)=0.0
    4 CONTINUE
      IF(FNUM.LT.0.)S(I+3)=180.-S(I+3)
      IF(SM(ID,ID).EQ.0.)E(ID)=0.
    1 CONTINUE
      IF(ABS(S(1)-S(2)).LT..00008)E(1)=E(2)
      DO 2 I=1,6
      SM(I,I)=E(I)
    2 V(I)=S(I)
      RETURN
      END
      SUBROUTINE STEEP(X,N,M,CHISQ,GRAD)
      include 'param.inc'
      DIMENSION Z(16),X(16),T(16),ZNT(16),SIGT(16),SIG(16)
      DIMENSION Y(16,16),XI(16),ZN(16)
      IF(N.EQ.0)RETURN
      CALL GRAD(X,Z)
      CALL CHISQ(X,G)
      DO 10 I=1,M
      DO 10 J=1,M
      Y(I,J)=0.0
      IF(I.EQ.J)Y(I,J)=1.0
   10 CONTINUE
      DO 200 K=1,N
      U=0
      DO 20 I=1,M
      T1=0
      DO 30 J=1,M
   30 T1=T1+Z(J)*Y(I,J)
      SIG(I)=T1
   20 U=U-SIG(I)*Z(I)
      ETA=AMIN1(1.0,ABS(2.*G/U))
      DO 40 I=1,M
   40 T(I)=X(I)-ETA*SIG(I)
      CALL GRAD(T,ZN)
      V=0
      DO 50 I=1,M
   50 V=V-ZN(I)*SIG(I)
      CALL CHISQ(T,G1)
      D=3.*ABS(G-G1)/ETA+U+V
      ARGU=D*D-U*V
      IF (ARGU.LT.0.0) ARGU=0.0
      W=SQRT(ARGU)
      ALPHA=ETA*(1.-(V+W-D)/(V-U+2.*W))
      DO 60 I=1,M
      SIG(I)=-1.0*ALPHA*SIG(I)
      X(I)=X(I)+SIG(I)
   60 ZN(I)=Z(I)
      CALL GRAD(X,Z)
      CALL CHISQ(X,G)
      T1=0
      DO 70 I=1,M
      XI(I)=Z(I)-ZN(I)
   70 T1=T1+XI(I)*SIG(I)
      DO 80 I=1,M
   80 SIGT(I)=SIG(I)/T1
      T1=0
      DO 90 I=1,M
      T2=0
      DO 100 J=1,M
  100 T2=T2+XI(J)*Y(I,J)
      ZN(I)=T2
   90 T1=T1+ZN(I)*XI(I)
      DO 110 I=1,M
      T2=0
      DO 120 J=1,M
  120 T2=T2+XI(J)*Y(J,I)
  110 ZNT(I)=T2/T1
      DO 131 I=1,M
      DO 130 J=1,M
  130 Y(I,J)=Y(I,J)+SIG(I)*SIGT(J)-ZN(I)*ZNT(J)
  131 CONTINUE
  200 CONTINUE
      RETURN
      END
      SUBROUTINE SUMMAT(IPM)
      include 'param.inc'
      REAL LAMDA
      INTEGER PTR
      COMMON/F1/RJAC(MSZ,MSZ),V1(MSZ)
      COMMON/REFLS/IREFS(IRS),REFS(IRS,8),FF(IRS),IC(NUPH)
      COMMON/G1/THMIN,STEP,THMAX,U,V,W,LAMDA(2),TMV,CTHM,RLIM,SB,WDT,
     *UU,VV,WW,SYNPOL
      COMMON/CNTRLS/JOBTYP,NPROF,NPHASE,IOT,IPL,IPC,MAT,NXT,LST1,LST2,
     *LST3,IPL1,IPL2,MCYCLE,EPS,MAXS,IFOUR,NFMT
      COMMON/PARAMS/GLB(16),AGLB(16),LGLB(16),PAR(NUPH,20),
     *APAR(NUPH,20),LPAR(NUPH,20),XL(NATS,11),LP(NATS,11),A(NATS,11),
     *PTR(NATS),RELAX(4),RATIO(2)
      CHARACTER*4 ATEXT(NATS),NTYP(NATS)
      COMMON/PARAC/ATEXT,NTYP
      COMMON/JNK/ALOW(30),AHIGH(30),POS(30),BCK(30),
     *NATOM(NUPH),PREF(NUPH,3),NBCKGD,NEXCRG,NSCAT,IWT
      CHARACTER SYMB(NUPH,20)*1
      CHARACTER TITLE*80,PHSNM(NUPH)*80
      COMMON/CHAR/SYMB,TITLE,PHSNM
      COMMON/COEFF/AC(10,16),POSI(30),SCAT(30),DFP(16),DFPP(16)
      CHARACTER*4 NAM(16)
      COMMON/COEFC/NAM
      COMMON/DATAX/Y(IDSZ),YC(IDSZ),KR(IDSZ),BK(IDSZ),NPTS,WT(IDSZ)
      COMMON/G3/COND,IORD1,IORD2,TH,NUM
      COMMON/G4/TANN(NIPOI),DERSTO(NIPOI,MSZ)
      COMMON/PRFX/TANTH,TL,GAM,PASYM,DERHW,DERGAM,DERASY,DERZ,DERTTB
      COMMON/RBRAGG/COBS(IRS),SOMEGA(IRS),VAR(IRS)
      DIMENSION DERIV(MSZ)
      DIMENSION DUMP(NIPOI),SOME(NIPOI)
      LOGICAL VERT
      DO 2 J=1,MAXS
    2 DERIV(J)=0.
      YCALC=0.
C-----CALCULATE THE CONTRIBUTION OF THE REFLECTIONS ORD1 TO ORD2 TO THE
C-----DERIVATIVES W.R.T. THE PROFILE INTENSITY YOBS
      IF(IORD1.EQ.0)GOTO 12
      DO 33 I=IORD1,IORD2
      IPH=IREFS(I)/256/256/256/4
      J=MOD(I,NIPOI)+1
      ICZ=MOD(IREFS(I)/256/256/256,4)
C-----VERT=.TRUE. IF RIETVELD ASYMMETRY CORRECTION IS TO BE APPLIED
      VERT=REFS(I,2).LE.RLIM.AND.NPROF.LT.4
      DELTA=TH-REFS(I,2)-GLB(1)
      DELT=DELTA*DELTA
      TL=REFS(I,1)
      GAM=REFS(I,4)
      TANTH=TANN(J)
      PASYM=PAR(IPH,14)
      IF(.NOT.VERT)GOTO 4
      YX=DELT*SIGN(1.,DELTA)
      Z=1.-PAR(IPH,14)*YX/TANN(J)
      IF(Z.LE.0.)Z=.0001
      GOTO 5
    4 Z=1.
    5 OMEGA=REFS(I,3)*Z*PROFIL(NPROF,DELTA)*PAR(IPH,1)
      IDUMP=I+1-IORD1
      DUMP(IDUMP)=OMEGA*FF(I)
      SOME(IDUMP)=OMEGA
      YCALC=YCALC+OMEGA*FF(I)
      IF(JOBTYP.GT.2)GOTO 33
      DO 3 K=1,MAXS
      DER=1.
      DO 6 M=3,5
      IF(LPAR(IPH,M).EQ.K) DER=DERHW
    6 CONTINUE
      X1=0.
      IF(VERT)X1=2.0*PAR(IPH,14)*SIGN(1.,DELTA)/TANN(J)/Z
      IF(LGLB(1).EQ.K) DER=DELTA*X1+DERZ
      IF(LPAR(IPH,14).EQ.K) DER=DERASY
      IF(LPAR(IPH,14).EQ.K.AND.VERT)DER=-YX/TANN(J)/Z
      IF(K.EQ.LPAR(IPH,15)) DER=DERGAM
      IF(K.EQ.LPAR(IPH,16)) DER=DERGAM
      IF(K.EQ.LPAR(IPH,17)) DER=DERGAM
    8 DO 7 M=6,11
      IF(LPAR(IPH,M).EQ.K) DER=DELTA*X1+DERTTB
    7 CONTINUE
    3 DERIV(K)=DERSTO(J,K)*DER*OMEGA+DERIV(K)
   33 CONTINUE
C-----FORM SUMS
   12 IF(NBCKGD.NE.0)GOTO 11
      DO 10 I=2,6
      IF(LGLB(I).EQ.0)GOTO 10
      KM=LGLB(I)
      IF(I.EQ.2.AND.TH.EQ.0.0) TH=1.0
      DERIV(KM)=DERIV(KM)+TH**(I-2)
   10 CONTINUE
      IF(LGLB(7).EQ.0) GO TO 101
      KM=LGLB(7)
      DERIV(KM)=DERIV(KM)+1.0/TH
  101 CONTINUE
   11 YC(IPM)=YCALC
      IF(IWT.EQ.1.AND.MCYCLE.GT.1)
     *  WT(IPM)=1.0/(YC(IPM)+BK(IPM))
      DELTA=Y(IPM)-BK(IPM)-YCALC
      DO 9 J=1,MAXS
      X=DERIV(J)*WT(IPM)
      V1(J)=V1(J)+DELTA*X
      DO 9 K=J,MAXS
   9  RJAC(J,K)=RJAC(J,K)+X*DERIV(K)
      IF(IORD1.EQ.0) GO TO 20
C     SET UP ARRAY OF CALCULATED OBSERVATIONS
      IF(JOBTYP.GT.2) WT(IPM)=1.0/(YC(IPM)+BK(IPM))
      DO 30 I=IORD1,IORD2
      IDUMP=I+1-IORD1
      SOMEGA(I)=SOMEGA(I)+SOME(IDUMP)
      IF(NBCKGD.NE.0) GO TO 29
      FRAC=DUMP(IDUMP)/(YCALC+BK(IPM))
      COBS(I)=COBS(I)+Y(IPM)*FRAC
      VAR(I)=VAR(I)+FRAC*FRAC/WT(IPM)
      GO TO 30
   29 IF(YCALC.EQ.0) GO TO 30
      FRAC=DUMP(IDUMP)/YCALC
      COBS(I)=COBS(I)+(Y(IPM)-BK(IPM))*FRAC
      VAR(I)=VAR(I)+FRAC*FRAC/WT(IPM)
   30 CONTINUE
   20 RETURN
      END
      SUBROUTINE SORT(IPHASE)
      include 'param.inc'
      COMMON/REFLS/IREFS(IRS),REFS(IRS,8),MLT(IRS),ICX(NUPH)
      DIMENSION L(IRS),TEMP(IRS),ITEMP(IRS)
      EQUIVALENCE (ITEMP(1),TEMP(1))
      IOF=0
      IF(IPHASE.GE.2) THEN
      DO 36 I=1,IPHASE-1
   36 IOF=IOF+ICX(I)
      ENDIF
      IC=ICX(IPHASE)
      DO 2 I=1,IC
    2 L(I)=I
      M=IC-1
      DO 15 LL=1,5000
      K=0
      DO 10 J=1,2
      DO 10 I=J,M,2
      INP=L(I)+IOF
      INP1=L(I+1)+IOF
      R=REFS(INP1,2)-REFS(INP,2)
      IF(R) 20,10,10
   20 IZ=L(I)
      L(I)=L(I+1)
      L(I+1)=IZ
      K=K+1
   10 CONTINUE
      IF(K.EQ.0) GO TO 5
   15 CONTINUE
    5 DO 30 I=1,IC
   30 ITEMP(I)=IREFS(I+IOF)
      DO 31 I=1,IC
      LX=L(I)
      IREFS(I+IOF)=ITEMP(LX)
   31 CONTINUE
      DO 33 J=1,4
      DO 32 I=1,IC
   32 TEMP(I)=REFS(I+IOF,J)
      DO 33 I=1,IC
      LX=L(I)
   33 REFS(I+IOF,J)=TEMP(LX)
      DO 34 I=1,IC
   34 ITEMP(I)=MLT(I+IOF)
      DO 35 I=1,IC
      LX=L(I)
   35 MLT(I+IOF)=ITEMP(LX)
      RETURN
      END
      SUBROUTINE LOOKUP(K,N,NSCAT,IXRAY,JOB)
      include 'param.inc'
      INTEGER PTR,TBXPTR
      CHARACTER*4 TABNC(128),TBXC(256)
      COMMON/TABCHR/TABNC,TBXC
      COMMON/TABLES/TBX(256,10),TBD(128,10),TABN(128)
      COMMON/COEFF/AC(10,16),POSI(30),SCAT(30),DFP(16),DFPP(16)
      CHARACTER*4 NAM(16)
      COMMON/COEFC/NAM
      COMMON/PARAMS/GLB(16),AGLB(16),LGLB(16),PAR(NUPH,20),
     *APAR(NUPH,20),LPAR(NUPH,20),XL(NATS,11),LP(NATS,11),A(NATS,11),
     *PTR(NATS),RELAX(4),RATIO(2)
      CHARACTER*4 ATEXT(NATS),NTYP(NATS)
      COMMON/PARAC/ATEXT,NTYP
      COMMON/DC/SAVE(NUPH,8),NSAVE
      COMMON/JNK/ALOW(30),AHIGH(30),POS(30),BCK(30),
     *NATOM(NUPH),PREF(NUPH,3),NBCKGD,NEXCRG,NSCATX,IWT
      CHARACTER SYMB(NUPH,20)*1
      CHARACTER TITLE*80,PHSNM(NUPH)*80
      COMMON/CHAR/SYMB,TITLE,PHSNM
      IOF=0
      NS=0
      IF(K.GT.1) THEN
      DO 121 I=1,K-1
  121 IOF=IOF+NATOM(I)
      ENDIF
      IF(K.GT.1)NS=NSAVE
      IF(JOB.EQ.2.OR.JOB.EQ.4) GOTO 70
      DO 10 I=1,N
      NS=MAX0(NSCAT,NS)
      DO 20 J=1,NS
      IF(NTYP(I+IOF).EQ.NAM(J)) GOTO 30
   20 CONTINUE
      DO 40 J=1,213
       IF(NTYP(I+IOF).EQ.TBXC(J)) GOTO 50
   40 CONTINUE
      WRITE(16,41)NTYP(I+IOF)
   41 FORMAT(1x,' SCATTERING COEFFICIENTS NOT FOUND FOR ',A4)
      STOP ' SCATTERING DATA MISSING'
   30 PTR(I+IOF)=J
      GOTO 10
   50 NS=NS+1
      PTR(I+IOF)=NS
      DO 60 L=1,9
   60 AC(L,NS)=TBX(J,L)
      NAM(NS)=TBXC(J)
      TBXPTR=TBX(J,10)+0.5
      DFP(NS)=TBD(TBXPTR,IXRAY)
      DFPP(NS)=TBD(TBXPTR,IXRAY+5)
   10 CONTINUE
      NSAVE=NS
      RETURN
   70 DO 80 I=1,N
      NS=MAX0(NSCAT,NS)
      DO 90 J=1,NS
      IF(NTYP(I+IOF).EQ.NAM(J))GOTO 100
   90 CONTINUE
      DO 110 J=1,86
      IF(NTYP(I+IOF).EQ.TABNC(J))GOTO 120
  110 CONTINUE
      WRITE(16,111)NTYP(I+IOF)
  111 FORMAT(34H SCATTERING LENGTHS NOT FOUND FOR ,A4)
      STOP ' SCATTERING DATA MISSING'
  100 PTR(I+IOF)=J
      GOTO 80
  120 NS=NS+1
      PTR(I+IOF)=NS
      DFP(NS)=TABN(J)
      NAM(NS)=TABNC(J)
   80 CONTINUE
      NSAVE=NS
      RETURN
      END
      SUBROUTINE CELL2
      include 'param.inc'
      COMMON/SPG/NSPGRP,K,FILL(13)
      COMMON/CELLX/A,B,C,ALPHA,BETA,GAMMA,AL(3,3)
      EQUIVALENCE (AL(1,1),ASTAR),(AL(2,2),BSTAR),(AL(3,3),CSTAR),
     *(AL(1,2),BACOSC),(AL(1,3),CACOSB),(AL(2,3),BCCOSA)
   35 FORMAT (22H THE LAUE SYMMETRY IS ,A8)
      CHARACTER*8 LAU(14)
      DATA LAU(1),LAU(2),LAU(3),LAU(4),LAU(5),LAU(6),LAU(7),LAU(8),
     *LAU(9),LAU(10),LAU(11),LAU(12),LAU(13),LAU(14)
     */'1BAR','2/M','MMM','4/M','4/MMM','3BAR   R','3BAR M R','3BAR',
     *'3BAR M 1','3BAR 1 M','6/M','6/MMM','M3','M3M'/
C     SIN1(A) = SIN(A*6.28318531)
      COS1(A) = COS(A*6.28318531)
      WRITE(16,35) LAU(NSPGRP)
      IF ( K.GT.3 ) STOP 5001
      GO TO (1002,1071,1072,1073,1073,1074,1074,1075,1075,1075,1075,
     1 1075,1076,1076),NSPGRP
 1071 IF ( K.NE.1 ) ALPHA=90.0
      IF ( K.NE.2 ) BETA=90.0
      IF ( K.NE.3) GAMMA=90.0
      GO TO 1002
 1073 B = A
 1072 ALPHA = 90.0
      BETA = 90.0
      GAMMA = 90.0
      GO TO 1002
 1074 B = A
      C = A
      BETA = ALPHA
      GAMMA = ALPHA
      GO TO 1002
 1075 B = A
      ALPHA = 90.
      BETA = 90.0
      GAMMA = 120.0
      GO TO 1002
 1076 B = A
      C = A
      GO TO 1072
 1002 COSA = COS1(ALPHA/360.0)
  105 COSB = COS1(BETA/360.0)
 1052 COSC = COS1(GAMMA/360.0)
 1053 SINA=SQRT(1.0-COSA**2)
      SINB=SQRT(1.0-COSB**2)
      SINC=SQRT(1.0-COSC**2)
      COSASR=(COSB*COSC-COSA)/(SINB*SINC)
      COSBSR=(COSA*COSC-COSB)/(SINA*SINC)
      COSCSR=(COSA*COSB-COSC)/(SINA*SINB)
      SINASR=SQRT(1.0-COSASR**2)
      SINBSR=SQRT(1.0-COSBSR**2)
      ASTAR=1.0/(A*SINC*SINBSR)
      BSTAR=1.0/(B *SINC*SINASR)
      CSTAR=1.0/(C*SINA*SINBSR)
      BCCOSA=BSTAR*CSTAR*COSASR*2.0
      CACOSB=CSTAR*ASTAR*COSBSR*2.0
      BACOSC=BSTAR*ASTAR*COSCSR*2.0
      ASTAR=ASTAR*ASTAR
      BSTAR=BSTAR*BSTAR
      CSTAR=CSTAR*CSTAR
      RETURN
      END
      SUBROUTINE PVOIT(X,Y,U,V,DUDX,DVDY)
      Y=1.201122409*X
      U=EXP(-X*X)
      DUDX=-2.0*X*U
      V=1.0/(1.0+Y*Y)
      DVDY=-2.0*Y*V*V
      RETURN
      END
      SUBROUTINE WERF1(X,Y,U,DUDX,DUDY)
      CALL FC01A(X,Y,U,V)
      IF (U.LT.1.0E-05) U=1.0E-05
      DUDX=-2.0*X*U + 2.0*Y*V
      DUDY=-1.128379167 + 2.0*X*V + 2.0*Y*U
      RETURN
      END
      SUBROUTINE FC01A(SSXSI,SSETA,SSU,SSV)
C     FROM HARWELL SUBROUTINE LIBRARY - SEE REPORT AERE-R9185
      IMPLICIT REAL*8(A-H,O-Z)
      REAL * 4 SSXSI,SSETA,SSU,SSV
      DIMENSION S(41),T(41),C(30),D(30)
      XSI=SSXSI
      ETA=SSETA
      IENT=1
      GO TO 77
      ENTRY FC01AD(DXSI,DETA,DU,DV)
      IENT=2
      XSI=DXSI
      ETA=DETA
   77 CONTINUE
      PI=3.1415926535898D0
      PJ=2.0D0/DSQRT(PI)
      X=DABS(XSI)
      Y=DABS(ETA)
      Z=DSQRT(X*X+Y*Y)
      IF(Z-0.9D0)1,2,2
    1 ASSIGN 100 TO I1
      GO TO 30
    2 IF(Z-1.1D0)3,3,6
    3 ASSIGN 4 TO I1
      GO TO 30
    4 U1 = U
      V1 = V
      ASSIGN 5 TO I1
      GO TO 40
    5 U2 = U
      V2 = V
      U=5.0D0*((Z-0.9D0)*U2+(1.1D0-Z)*U1)
       V=5.0D0*((Z-0.9D0)*V2+(1.1D0-Z)*V1)
      GO TO 100
    6 IF(Z-3.9D0)7,120,120
    7 ASSIGN 99 TO I1
      ASSIGN 100 TO I2
    8 IF(Y-1.3D0)40,9,9
    9 IF(Y-1.5D0)10,10,13
   10 ASSIGN 11 TO I1
      GO TO 40
   11 U1 = U
      V1 = V
      ASSIGN 12 TO I1
       H=0.3D0
      GO TO 54
   12 U2 = U
      V2 = V
      U=5.0D0*((Y-1.3D0)*U2+(1.5D0-Y)*U1)
      V=5.0D0*((Y-1.3D0)*V2+(1.5D0-Y)*V1)
      GO TO 99
   13 IF(Y-1.8D0)14,15,15
   14 H=0.3D0
      GO TO 54
   15 IF(Y-2.0D0)16,16,101
   16 ASSIGN 17 TO I1
      H=0.3D0
      GO TO 54
   17 ASSIGN 18 TO I1
      U1 = U
      V1 = V
        H=0.4D0
      GO TO 54
   18 U2 = U
      V2 = V
      U=5.0D0*((Y-1.8D0)*U2+(2.0D0-Y)*U1)
        V=5.0D0*((Y-1.8D0)*V2+(2.0D0-Y)*V1)
      GO TO 99
  101 IF(Y-2.5D0)102,103,103
  102 ASSIGN 99 TO I1
      H=0.4D0
      GO TO 54
  103 IF(Y-2.7D0)104,104,107
  104 ASSIGN 105 TO I1
      H=0.4D0
      GO TO 54
  105 U1 = U
      V1 = V
      ASSIGN 106 TO I1
      H=0.5D0
      GO TO 54
  106 U2 = U
      V2 = V
      U=5.0D0*((Y-2.5D0)*U2+(2.7D0-Y)*U1)
      V=5.0D0*((Y-2.5D0)*V2+(2.7D0-Y)*V1)
      GO TO 99
  107 IF(Y-3.75D0)108,109,109
  108 ASSIGN 99 TO I1
      H=0.5D0
      GO TO 54
  109 IF(Y-3.85D0)112,113,113
  112 ASSIGN 110 TO I1
      H=0.5D0
      GO TO 54
  110 U1 = U
      V1 = V
      ASSIGN 111 TO I1
      H=0.6D0
      GO TO 54
  111 U2 = U
      V2 = V
       U=10.0D0*((Y-3.75D0)*U2+(3.85D0-Y)*U1)
      V=10.0D0*((Y-3.75D0)*V2+(3.85D0-Y)*V1)
      GO TO 99
  113 ASSIGN 99 TO I1
      H=0.6D0
      GO TO 54
  120 IF(Z-4.1D0)121,121,123
  121 ASSIGN 99 TO  I1
      ASSIGN 122 TO I2
      GO TO 8
  122 U1 = U
      V1 = V
      ASSIGN 130 TO I2
      GO TO 125
  123 IF(Z-5.8D0)124,124,131
  124 ASSIGN 99 TO I1
      ASSIGN 100 TO I2
  125 A2=2.0D0*X*Y
      A6 = X*X+Y*Y
      GO TO 64
  130 U2 = U
      V2 = V
      U=5.0D0*((4.1D0-Z)*U1+(Z-3.9D0)*U2)
      V=5.0D0*((4.1D0-Z)*V1+(Z-3.9D0)*V2)
      GO TO 100
  131 IF(Z-6.2D0)134,134,132
  132 ASSIGN 100 TO I2
  133 A2=2.0D0*X*Y
      A6 = X*X+Y*Y
      ASSIGN 99 TO I1
      GO TO 61
  134 ASSIGN 135 TO I2
      GO TO 133
  135 U1 = U
      V1 = V
      ASSIGN 136 TO I2
      GO TO 125
  136 U2 = U
      V2 = V
      U=2.5D0*((Z-5.8D0)*U1+(6.2D0-Z)*U2)
      V=2.5D0*((Z-5.8D0)*V1+(6.2D0-Z)*V2)
      GO TO 100
   30 S(1)=X
      T(1)=Y
      S(2)=X*X-Y*Y
      T(2)=2.0D0*X*Y
      E1=DEXP(-S(2))
      E2=DCOS(T(2))
      E3=DSIN(T(2))
      J=1
      U=1.0D0-PJ*T(1)
      V=PJ*S(1)
      E4=E1*(E2*U+E3*V)
      V=E1*(E2*V-E3*U)
      U=E4
      DO 33 I=3,41,2
      S(I)=(S(I-2)*S(2)-T(I-2)*T(2))/J
      T(I)=(S(I-2)*T(2)+T(I-2)*S(2))/J
      J=J+1
      UE=-PJ*T(I)/I
      VE=PJ*S(I)/I
      E5=E1*(E2*UE+E3*VE)
      VE=E1*(E2*VE-E3*UE)
      UE=E5
      IF(DABS(UE)-.1D-6)31,31,32
   31 IF(DABS(VE)-.1D-6)34,34,32
   32 U=U+UE
      V=V+VE
   33 CONTINUE
   34 GO TO 80
   40 S(1)=0.46752093659378D0
      T(1)=0.17912253437764D0
      S(2)=0.20798989547526D0
      T(2)=0.29084679383680D0
      S(3)=0.07508765293733D0
      T(3)=0.21954582524901D0
      S(4)=0.03526344064839D0
      T(4)=0.16054790355559D0
      I=AMIN0(3,IDINT(X))
      A=I+0.5
      B=0.7D0
      C(1)=S(I+1)
      D(1)=T(I+1)
      C(2)=-2.0D0*A*C(1)+2.0D0*B*D(1)
      D(2)=-2.0D0*A*D(1)-2.0D0*B*C(1)+PJ
      F=X-A
      G=Y-B
      U=C(1)+C(2)*F-D(2)*G
      V=D(1)+C(2)*G+D(2)*F
      FA=F*F-G*G
      GA=2.0D0*F*G
      DO 43 I=3,30
      C(I)=(-2.0D0*A*C(I-1)+2.0D0*B*D(I-1)-2.0D0*C(I-2))/(I-1)
      D(I)=(-2.0D0*A*D(I-1)-2.0D0*B*C(I-1)-2.0D0*D(I-2))/(I-1)
      UB=C(I)*FA-D(I)*GA
      VB=C(I)*GA+D(I)*FA
      FB=FA*F-GA*G
      GA=FA*G+GA*F
      FA=FB
      IF(DABS(UB)-.1D-6)41,41,42
   41 IF(DABS(VB)-.1D-6)44,44,42
   42 U=U+UB
      V=V+VB
   43 CONTINUE
   44 GO TO 80
   54 E2=DEXP(H*H)
      E3=1.0D0/(E2*E2)
      A1=PI*(X*X+Y*Y)
      U=Y*H/A1
      V=X*H/A1
      E1=1.0D0
      W=0.0D0
      DO 56 L=1,50
      W=W+H
      E1=E1*E2*E3**L
      A1=X-W
      A2=X+W
      A3=(A1*A1+Y*Y)*PI
      A4=(A2*A2+Y*Y)*PI
      A5=H*E1*Y*(1.0D0/A3+1.0D0/A4)
      A6=H*E1*(A1/A3+A2/A4)
      U=U+A5
      V=V+A6
      IF(DABS(A5)-.1D-6)55,55,56
   55 IF(DABS(A6)-.1D-6)57,57,56
   56 CONTINUE
   57 GO TO 80
   61 ALP1=0.01995324205905D0/PI
      X1=2.020182870456086D0
      I=1
      GO TO 72
   62 U=UA
      V=VA
      ALP1=0.3936193231522D0/PI
      X1=0.958572464613819D0
      I=2
      GO TO 72
   63 ALP2=0.9453087204829D0/PI
      U=U+UA+Y*ALP2/A6
      V=V+VA+X*ALP2/A6
      GO TO 80
   64 ALP1=.7640432855233D-5/PI
      X1=3.436159118837738D0
      I=3
      GO TO 72
   65 U=UA
      V=VA
      X1=2.532731674232790D0
      ALP1=.1343645746781D-2/PI
      I=4
      GO TO 72
   66 U=U+UA
      V=V+VA
      ALP1=.3387439445548D-1/PI
      X1=1.756683649299882D0
      I=5
      GO TO 72
   67 U=U+UA
      V=V+VA
      ALP1=0.2401386110823D0/PI
      X1=1.036610829789514D0
      I=6
      GO TO 72
   68 U=U+UA
      V=V+VA
      ALP1=0.6108626337353D0/PI
      X1=0.342901327223705D0
      I=7
      GO TO 72
   69 U=U+UA
      V=V+VA
      GO TO 80
   72 A1=X*X-Y*Y-X1*X1
      A3=A1*A1+A2*A2
      A4=X*A2-Y*A1
      A5=X*A1+Y*A2
      UA=2.0D0*ALP1*A4/A3
      VA=2.0D0*ALP1*A5/A3
      GO TO (62,63,65,66,67,68,69),I
   80 IF(X-XSI)84,81,84
   81 IF(Y-ETA)82,90,82
   82 E1=DEXP(Y*Y-X*X)
      E2=2.0D0*X*Y
      E3=DCOS(E2)
      U=-U+2.0D0*E1*E3
      E3=DSIN(E2)
      V=-V-2.0D0*E1*E3
      IF(X-XSI)90,83,90
   83 V=-V
      GO TO 90
   84 IF(Y-ETA)82,83,82
   90 GO TO I1,(4,5,11,12,17,18,99,100,105,106,110,111)
   99 GOTOI2,(100,122,130,135,136)
  100 IF(IENT.NE.1)GO TO 220
      SSU=U
      SSV=V
      GO TO 222
  220 DU=U
      DV=V
  222 RETURN
      END
