C     Last change:  AMF  29 Mar 95    1:36 pm
C
C
C     PROGRAM FORDAP(INPUT,OUTPUT,PUNCH,TAPE60=INPUT,TAPE61=OUTPUT,
C    1TAPE62=PUNCH,TAPE1,TAPE2,TAPE8)
C     ZALKIN FORDAP PROGRAM MODIFIED BY IBERS FOR BNL SYSTEM AND
C     ASSOCIATED PROGRAMS , AUGUST,1962.
C     VIA NORTHWESTERN UNIVERSITY AND UNIVERSITY OF CALIFORNIA,IRVINE.
C
C     MODIFIED FOR 16K 360/44,WITH TWIN 2315 DISK DRIVES,AT THE DEPT.
C     OF CHEMISTRY,UNIVERSITY OF CANTERBURY,CHRISTCHURCH,NEW ZEALAND.
C     VERSION EXPORTED OCTOBER 1968
C


      CHARACTER*20 NAME2,NAME60,NAME61
 2000 format (a,$)
 2001 format (a)



      DIMENSION TRGTBL(1501),ISST(10)
      DIMENSION NDATA( 5900),WORDS(24)
      DIMENSION DUMMY( 6486),U(3,3)
      COMMON TRGTBL,NORIEN,NXX,NYX,NZ,XMIN,XMAX,YMIN,YMAX,ZMIN
      COMMON ZMAX,DX,DY,DZ,AMAX,BMAX,JMAX,KMAX,KMIN,NMEMMX,NEGINA
      COMMON NCD,NEF,NEFTBL,NX,NY,NCUE,N1,RMAX,NLAST,NFPLT,YFL
      COMMON N11,IN,III,IVXX,IXX,NN2,NN3,NN5,ISST,A,B,C,COSA
      COMMON COSB,COSC,NSEC,NGRID,NSKIP,CNTRNO,CNTNEG,NEXPAN
      COMMON WORDS,NDATA
      COMMON/PEAK/RMIN,MWG,WGM,WGN,WGO
      EQUIVALENCE (DUMMY(1),TRGTBL(1))


      write (*, 2000) 'INPUT DATA IS ON FILE ?.fou '
      read (*, 2001) NAME2
      I = index (NAME2, '.')
      if (I .eq. 0) I = index (NAME2, ' ')
      NAME2(I:I) = '.'
      NAME2 = NAME2(:I) // 'fou'

      write (*,2000) 'INPUT PARAMETERS ARE ON FILE ?.inp '
      read (*,2001) NAME60
      I = index (NAME60, '.')
      if (I .eq. 0) I = index (NAME60, ' ')
      NAME60(I:I) = '.'
      NAME60 = NAME60(:I) // 'inp'

      write (*,2000) 'OUTPUT PARAMETERS GO TO FILE ?.map '
      read (*,2001) NAME61
      I = index (NAME61, '.')
      if (I .eq. 0) I = index (NAME61, ' ')
      NAME61(I:I) = '.'
      NAME61 = NAME61(:I) // 'map'

      OPEN (UNIT=2,FILE=NAME2,ERR=99,STATUS='OLD',FORM='UNFORMATTED')
      OPEN (UNIT=60,FILE=NAME60,ERR=98,STATUS='OLD')
      OPEN (UNIT=61,FILE=NAME61,ERR=97,STATUS='UNKNOWN')
      OPEN (UNIT=16,ERR=96,STATUS='SCRATCH')
      OPEN (UNIT=1,ERR=96,STATUS='SCRATCH')
      OPEN (UNIT=8,ERR=96,STATUS='SCRATCH')

      GO TO 210

   99 STOP 'ERROR OPENING UNIT 2'
   98 STOP 'ERROR OPENING UNIT 60'
   97 STOP 'ERROR OPENING UNIT 61'
   96 STOP 'ERROR OPENING SCRATCH FILE'

  110 IF(ISST(1).EQ.0) GO TO 210
      REWIND IVXX
  210 DO 1 I=1,6486
 1    DUMMY(I)=0.0
*$ifdef debug
      write( *,* ) 'Opened files' 
*$endif      
C
C     READ PARAMETER CARDS DEFINING CALCULATION AND DATA FOR FOURIER
      CALL SETUP(VOL)
*$ifdef debug
      write( *,* ) 'Returned from setup' 
*$endif      
      
      IF(ISST(2))63,65,63
C     TRANSFORM UNIT CELL AND INDICES FOR FOURIER IN A GENERAL PLANE
   63 CALL TRANS(U,DETU)
*$ifdef debug
      write( *,* ) 'Returned from trans'
*$endif      
 65   IF (ISST(1))42,42,64
C     EXPAND AND PACK THE REFLECTION DATA FOR A FOURIER SUMMATION
   64 CALL RIP(U,VOL)
*$ifdef debug
      write( *,* ) 'Returned from rip' 
*$endif      
C     CARRY OUT THE SPECIFIED FOURIER SUMMATION
      CALL SIG100
*$ifdef debug
      write( *,* ) 'Returned from sig'
*$endif      
   42 IF (ISST(3))43,44,43
C     PRINT THE NUMERIC ANSWERS FOR THE FOURIER SUMMATION
   43 CALL FRPRNT
*$ifdef debug
      write( *,* ) 'Returned from frprnt'
*$endif      
   44 IF(ISST(5).EQ.0) GO TO 46
C     PLOT AN ALPHAMERIC FOURIER
C
   47 CALL FRPLOT
*$ifdef debug
      write( *,* ) 'Returned from frplot' 
*$endif      
   46 IF(ISST(4).EQ.0) GO TO 110
C     SORT AND LIST THE MAXIMUM POINTS IN ORDER OF MAGNITUDE
   45 CALL FRPEAK(U,DETU,VOL)
*$ifdef debug
      write( *,* ) 'Returned from frpeak'
*$endif      
   
      CALL IPEAK(U)
*$ifdef debug
      write( *,* ) 'Returned from ipeak'
*$endif      
      
      END



      SUBROUTINE CHEKUP(XMIN,XMAX,NORIEN,N,DX)
    1 DIFFX=XMAX-XMIN
      IF (DIFFX)3,20,5
    3 XMIN=XMAX
      XMAX=XMIN-DIFFX
      GO TO 1
    5 IF (DIFFX-1.)30,30,10
   10 XMAX=.5
      XMIN=-.5
      GO TO 1
   50 DIFFX=AINT(DIFFX/DX+0.001)*DX
      XMAX=XMIN+DIFFX
      GO TO 60
   20 NORIEN=N
      DX=.001
      GO TO 60
   30 IF (XMAX-1.)50,50,40
   40 XMIN=-.5
      XMAX=XMIN+DIFFX
      GO TO 1
   60 RETURN
      END


      SUBROUTINE CROSS(BB,A,B,C)
      DIMENSION A(3),B(3)
      DIMENSION C(3),AI(6)
      DIMENSION BI(6),CI(3)
      DIMENSION BB(3,3)
      DO 115 I=1,3
      AI(I)=A(I)
      AI(I+3)=A(I)
      BI(I)=B(I)
  115 BI(I+3)=B(I)
      DO 119 I=1,3
  119 CI(I)=AI(I+1)*BI(I+2)-AI(I+2)*BI(I+1)
      CALL MV(BB,CI,C)
      RETURN
      END


      SUBROUTINE FRPEAK(U,DETU,VOL)
      INTEGER   BORDER
      DIMENSION BORDER(2)
      REAL*8 PEAKL,PEAKL1
      DIMENSION PEAKL(2),PEAKL1(2)
C     DIMENSION INEG(2)
      DIMENSION XX( 704),YY( 704)
      DIMENSION NXPREV(250),RRPREV(250)
      DIMENSION RR( 704),TRGTBL(1514)
      DIMENSION ISST(10),RRR(250)
      DIMENSION ZZ( 704),NNX( 704)
      DIMENSION NNY( 704),NNZ( 704)
      DIMENSION NNNX(250),NNNY(250)
      DIMENSION NYPREV(250),EXCESS(12)
      DIMENSION WORDS(24),SPACE(5)
      DIMENSION XKIM(3),R( 105,3)
      DIMENSION U(3,3)
      INTEGER   NNX,NNY,NNZ,NNNX,NNNY,NXPREV,NYPREV
      COMMON TRGTBL,AMAX,BMAX,EXCESS,RMAX,NLAST,YFPLT,YFL,N11
      COMMON IN,III,IVXX,IXX,NN2,NN3,NN5,ISST,A,B,C,COSA,COSB
      COMMON COSC,SPACE,NEXPAN,WORDS,R,XX,YY,ZZ,NNX,NNY,NNZ,RR
      COMMON/PEAK/ RMIN,NPEAK
      EQUIVALENCE (TRGTBL(1),RRPREV(1)),(TRGTBL(251),RRR(1)),
     $ (TRGTBL(501),NXPREV(1)),(TRGTBL(751),NYPREV(1)),
     $ (TRGTBL(1001),NNNX(1)), (TRGTBL(1251),NNNY(1))
       DATA BORDER(1),BORDER(2) /'  ',' +'/
       DATA IWRR /'00  '/
       DATA PEAKL(1),PEAKL(2) /'POSITIVE','NEGATIVE'/
       DATA PEAKL1(1),PEAKL1(2) /' PEAK','TROUGH'/
 1001 FORMAT (  '0 FRPEAK DOES NOT WORK ON SINGLE SECTIONS')
  406 FORMAT (  '0FOURPEAK AS WRITTEN FILTERS A MAXIMUM OF 250 PEAKS PER
     1 SECTION. (SECTION'I4' )')
 516  FORMAT (46H0FRPEAK IS LIMITED TO A SEARCH FOR  700 PEAKS.)
 703  FORMAT (///'1LIST OF '2A8    'POSITIONS AND THEIR HEIGHTS (THOSE O
     AN BORDER ARE MARKED +)'//)
1703  FORMAT('ED=ELECTRON DENSITY (INE/ACUBED IF INPUT DATA ARE IN'
     1'ELECTRONS)=(2./VOL)X(RMAX/999.)XRHO'/ ' IF LISTED, X(A) ETC REFER
     1 TO NEW CELL B WITH B2 PARALLEL A2' ,     ' , B1 PARALLEL TO A2XA3
     1 , B3 PARALLEL TO B1XB2'//)
  716 FORMAT('   NO.'9X'X'9X'Y'9X'Z'8X'RHO'5X'ED(E/AEX3)')
 7193 FORMAT (1H+,73X,24H    OLDX    OLDY    OLDZ)
 9106 FORMAT (1H+,73X,'    X(A)    Y(A)    Z(A)        D(000)')
 9107 FORMAT (1H0)
  810 FORMAT('0'I4,4X3F10.4,I8,3X'('F12.3')'A2)
 8107 FORMAT (1H+73X3F8.3,5XF8.3)
      INEG=0
      NPEAK=0
      REWIND IXX
      RRRMAX=RMAX
      REWIND IVXX
      RMIN=RMIN*RMAX/999.0
   10 LP=0
      JLAST=0
      IF(INEG.EQ.0) GO TO 4
      LT=4*ISST(8)
      GO TO 53
    4 LT=40*ISST(4)
   53 L=0
      NZ=0
      DO 70 M=1,250
      RRR(M)=RMIN
      NNNX(M)=0
      NNNY(M)=0
      RRPREV(M)=RMIN
      NXPREV(M)=0
   70 NYPREV(M)=0
      NUMM=1
      ASSIGN 100 TO NEWDAT
      ASSIGN 250 TO INDATA
 2110 READ(IVXX)CODE,RMAX,
     1           NX,NY,XMIN,XMAX,DELTAX,ZMAX,Z,NXIS,YMIN,YMAX,DY,AMAX,
     1BMAX,DZ,ZMIN,NZNZNZ
      NUM=0
      IF(CODE) 2111,200,2009
 2111 ASSIGN 700 TO NEWDAT
 2112 ASSIGN 240 TO INDATA
      GO TO 200
 2009 IF(NUMM.EQ.1) GO TO 11150
      ASSIGN 640 TO NEWDAT
      GO TO 2112
  100 NUM=NUM+1
      IF(NUM-NY) 200,200,2110
11150 NUMM=2
      IF(ZMAX.NE.ZMIN) GO TO 2
 1000 WRITE (III,1001)
      RETURN
 2    NXX=NX+2
      DO 120 K=1,3
      DO 120 I=1,NXX
  120 R(I,K)=RMIN
      NZ=NZ+1
      NXXX=NX+1
      J=0
      NYN=0
      JOF=0
      GO TO 100
  200 DO 220 I=1,NXX
      R(I,1)=R(I,2)
  220 R(I,2)=R(I,3)
      Y=Y+DY
      GO TO INDATA,(250,240)
  240 DO 245 I=2,NXXX
  245 R(I,3)=RMIN
      Z=Z-DZ
      ASSIGN 250 TO INDATA
      IF(JOF) 591,290,591
  250 READ(IVXX) CODE,CODE,
     1           Y,(R(I,3),I=2,NXXX)
      IF(JOF.NE.0) GO TO 591
      IF(INEG.EQ.0) GO TO 253
      DO 252 I=2,NXXX
  252 R(I,3)=-R(I,3)
  253 Y=Y-DY
  290 X=XMIN-DELTAX
      NYN=NYN+1
  305 DO 590 I=2,NXXX
      X=X+DELTAX
      IF(R(I,2).LE.RMIN) GO TO 590
      JJ=I-1
      JJJ=I+1
      DO 320 M=JJ,JJJ
      DO 320 K=1,3
      IF (R(I,2)-R(M,K))590,320,320
  320 CONTINUE
  400 J=J+1
      IF (J-250)408,408,405
  405 WRITE (III,406) NZ
      JOF=1
      J=J-1
      GO TO 591
  408 RRR(J)=R(I,2)
      NNNX(J)=I-1
      NNNY(J)=NYN
      IF(JLAST.EQ.0) GO TO 440
  409 DO 430 M=1,JLAST
      NEWY=IABS(NYN-NYPREV(M))
      IF (NEWY-1)410,410,430
  410 NEWX=IABS(I-1-NXPREV(M))
      IF (NEWX-1)415,415,430
  415 IF (R(I,2)-RRPREV(M))590,440,440
  430 CONTINUE
  440 IF(L.EQ.0) GO TO 510
      DO 480 K=1,L
      IF (NZ-NNZ(K)-1)480,450,480
  450 NEWY=IABS(NYN-NNY(K))
      IF (NEWY-1)460,460,480
  460 NEWX=IABS(I-1-NNX(K))
      IF (NEWX-1)470,470,480
  470 IF (R(I,2)-RR(K))590,475,475
  475 XX(K)=X
      YY(K)=Y
      ZZ(K)=Z
      RR(K)=R(I,2)
      NNX(K)=I-1
      NNY(K)=NYN
      NNZ(K)=NZ
      GO TO 590
  480 CONTINUE
  510 L=L+1
      IF (L- 700)520,520,515
  515 WRITE (III,516)
      WRITE(III,92537) NZ
92537 FORMAT('0+++PRESENT SECTION IS NUMBER',I4//)
      IF(RRRMAX.EQ.0.0) GO TO 700
      RMAX=RRRMAX
      GO TO 700
  520 XX(L)=X
      YY(L)=Y
      ZZ(L)=Z
      RR(L)=R(I,2)
      NNX(L)=I-1
      NNY(L)=NYN
      NNZ(L)=NZ
  590 CONTINUE
  591 GO TO NEWDAT,(100,640,700)
  640 ASSIGN 100 TO NEWDAT
      JLAST=J
      Z=Z+DZ
      DO 680 I=1,JLAST
      RRPREV(I)=RRR(I)
      NXPREV(I)=NNNX(I)
  680 NYPREV(I)=NNNY(I)
      GO TO 11150
  700 REWIND IVXX
      WRITE (III,703) PEAKL(INEG+1) , PEAKL1(INEG+1)
      IF(INEG.EQ.0) WRITE(III,1703)
  705 WRITE (III,716)
  719 IF (ISST(2))7192,7191,7192
 7192 WRITE (III,7193)
      GO TO 7194
 7191 DETU=1.0
      WRITE (III,9106)
      CALL SETU(U,A,B,C,COSA,COSB,COSC)
 7194 CONST=999./RMAX
17725 DO 740 K=1,L
      IF (RR(K)-RMIN)740,740,730
  730 ICE=K
      RBIG=RR(K)
      GO TO 760
  740 CONTINUE
      GO TO 809
  760 DO 780 K=1,L
      IF (RBIG-RR(K))770,780,780
  770 RBIG=RR(K)
      ICE=K
  780 CONTINUE
      RR(ICE)=RMIN
      NRR=RBIG*CONST+.5
      RBED=2.0*RBIG*DETU/VOL
 30   IF (ABS(XX(ICE)-XMAX)-0.0005)725,725,911
 911  IF (ABS(XX(ICE)-XMIN)-0.0005)725,725,912
 912  IF (ABS(YY(ICE)-YMAX)-0.0005)725,725,913
 913  IF (ABS(YY(ICE)-YMIN)-0.0005)725,725,914
 914  IF (ABS(ZZ(ICE)-ZMAX)-0.0005)725,725,915
 915  IF (ABS(ZZ(ICE)-ZMIN)-0.0005)725,725,916
  725 IBOR=2
      GO TO 1916
  916 IBOR=1
 1916 LP=LP+1
      IF(INEG.EQ.0) GO TO 922
      NRR=-NRR
      RBED=-RBED
  922 IF (NXIS-10)8103,8104,8105
 8104 CALL ROTATE(YY(ICE),ZZ(ICE))
      GO TO 8103
 8105 CALL ROTATE(XX(ICE),ZZ(ICE))
 8103 WRITE (III,810)LP,XX(ICE),YY(ICE),ZZ(ICE),NRR,RBED,BORDER(IBOR)
C     LL=IWRR+MOD(LP,10)*68719476736
C    A   +(MOD(LP,100)/10)*4398046511104
      LL=IWRR+MOD(LP,10)*65536 + (MOD(LP,100)/10)*16777216
      WRITE(IXX) WORDS(1), LL ,  XX(ICE),YY(ICE),ZZ(ICE)
     A  , IBOR,NNX(ICE),NNY(ICE),NNZ(ICE)
      NPEAK=NPEAK+1
      DO 8106 KIM1=1,3
 8106 XKIM(KIM1)=U(1,KIM1)*XX(ICE)+U(2,KIM1)*YY(ICE)+U(3,KIM1)*ZZ(ICE)
      IF (ISST(2))9108,9109,9108
 9109 DPAT=SQRT(XKIM(1)**2+XKIM(2)**2+XKIM(3)**2)
      GO TO 9110
 9108 DPAT=0
 9110 WRITE (III,8107)(XKIM(KIM1),KIM1=1,3),DPAT
 8102 IF (LP-LT)17725,809,809
  809 IF(INEG.GT.0) GO TO 880
      IF(ISST(8).EQ.0) GO TO 880
      REWIND IVXX
      INEG=1
      GO TO 10
 880  RETURN
      END


      SUBROUTINE EDINT(X,Y,Z,PSV,XX,IBOR)
      DIMENSION XX(3),PSV(3,3,3),P(20),X1(20),Y1(20),Z1(20)
      REAL A,B, C
      DIMENSION A(20,20),B(3),C(11,11)
      COMMON DUMMY(4440),A,B,C,P,X1,Y1,Z1
      COMMON /PEAK/ RMIN,NPEAK,DX,DY,DZ
C      ELECTRON DENSITY INTERPOLATION    ROLLETT PAGE 35
      X1( 1)=X
               Y1( 1)=Y+DY
                        Z1( 1)=Z-DZ
      X1( 2)=X-DX
               Y1( 2)=Y
                        Z1( 2)=Z-DZ
      X1( 3)=X
               Y1( 3)=Y
                        Z1( 3)=Z-DZ
      X1( 4)=X+DX
               Y1( 4)=Y
                        Z1( 4)=Z-DZ
      X1( 5)=X
               Y1( 5)=Y-DY
                        Z1( 5)=Z-DZ
      X1( 6)=X-DX
               Y1( 6)=Y+DY
                        Z1( 6)=Z
      X1( 7)=X
               Y1( 7)=Y+DY
                        Z1( 7)=Z
      X1( 8)=X+DX
               Y1( 8)=Y+DY
                        Z1( 8)=Z
      X1( 9)=X-DX
               Y1( 9)=Y
                        Z1( 9)=Z
      X1(10)=X
               Y1(10)=Y
                        Z1(10)=Z
      X1(11)=X+DX
               Y1(11)=Y
                        Z1(11)=Z
      X1(12)=X-DX
               Y1(12)=Y-DY
                        Z1(12)=Z
      X1(13)=X
               Y1(13)=Y-DY
                        Z1(13)=Z
      X1(14)=X+DX
               Y1(14)=Y-DY
                        Z1(14)=Z
      X1(15)=X
               Y1(15)=Y+DY
                        Z1(15)=Z+DZ
      X1(16)=X-DX
               Y1(16)=Y
                        Z1(16)=Z+DZ
      X1(17)=X
               Y1(17)=Y
                        Z1(17)=Z+DZ
      X1(18)=X+DX
               Y1(18)=Y
                        Z1(18)=Z+DZ
      X1(19)=X
               Y1(19)=Y-DY
                        Z1(19)=Z+DZ
C
      P( 1)=PSV(2,3,1)
      P( 5)=PSV(2,1,1)
      P( 2)=PSV(1,2,1)
      P( 3)=PSV(2,2,1)
      P( 4)=PSV(3,2,1)
      P( 6)=PSV(1,3,2)
      P( 7)=PSV(2,3,2)
      P( 8)=PSV(3,3,2)
      P( 9)=PSV(1,2,2)
      P(10)=PSV(2,2,2)
      P(11)=PSV(3,2,2)
      P(12)=PSV(1,1,2)
      P(13)=PSV(2,1,2)
      P(14)=PSV(3,1,2)
      P(15)=PSV(2,3,3)
      P(16)=PSV(1,2,3)
      P(17)=PSV(2,2,3)
      P(18)=PSV(3,2,3)
      P(19)=PSV(2,1,3)
C
      PSML=P(1)
      DO 30 I=2,19
   30 IF(PSML.GT.P(I)) PSML=P(I)
      IF(PSML.GT.0.0)     GO TO 32
      PSML=ABS(PSML) + 1.0
      DO 33 I=1,19
   33 P(I)=P(I)+PSML
C
C     SET UP LEAST-SQUARES MATRIX
C
   32 DO 4 I=1,19
      A(I, 1)=1.D0
      A(I, 2)=X1(I)
      A(I, 3)=Y1(I)
      A(I, 4)=Z1(I)
      A(I, 5)=X1(I)*X1(I)
      A(I, 6)=Y1(I)*Y1(I)
      A(I, 7)=Z1(I)*Z1(I)
      A(I, 8)=Y1(I)*Z1(I)
      A(I, 9)=X1(I)*Z1(I)
    4 A(I,10)=X1(I)*Y1(I)
      DO 5 I=1,11
      DO 5 J=1,11
    5 C(I,J)=0.D0
      DO 6 I=1,10
      DO 6 J=1,10
      DO 6 K=1,19
    6 C(I,J)=C(I,J)+A(K,I)*A(K,J)
      DO 7 I=1,10
      DO 7 J=1,19
    7 C(I,11)=C(I,11)+A(J,I)*ALOG(P(J))
C
C    INVERT AND SOLVE
C
      DO 8 LD=1,10
      C(11,11)=1.D0/C(1,1)
      DO 9 I=1,10
    9 C(11,I)=C(1,I+1)*C(11,11)
      DO 10 I=1,9
      DO 11 J=1,11
   11 C(I,J)=-C(I+1,1)*C(11,J)
      DO 10 J=1,10
   10 C(I,J)=C(I,J)+C(I+1,J+1)
      DO 8 J=1,11
    8 C(10,J)=C(11,J)
C   THE COEFFICIENTS ARE C(I,1) I=1,10   USE THEM TO SOLVE FOR X,Y,Z
C
      A(1,1)=-2.D0*C(5,1)
      A(1,2)=-C(10,1)
      A(1,3)=-C(9,1)
      A(2,1)=-C(10,1)
      A(2,2)=-2.D0*C(6,1)
      A(2,3)=-C(8,1)
      A(3,1)=-C(9,1)
      A(3,2)=-C(8,1)
      A(3,3)=-2.D0*C(7,1)
      B(1)=C(2,1)
      B(2)=C(3,1)
      B(3)=C(4,1)
      CALL GEMI(A,B,XX,3,NER)
      IF(ABS(XX(1)-X).GT.DX) GO TO 100
      IF(ABS(XX(2)-Y).GT.DY) GO TO 100
      IF(ABS(XX(3)-Z).GT.DZ) GO TO 100
      GO TO 110
  100 IBOR=3
  110 RETURN
      END


      SUBROUTINE FRPLOT
      DIMENSION ALPHA1(2,12),ALPHA2(2,12),ALPHA3(2,12)
      DIMENSION X( 103),ISST(10)
      DIMENSION Y(101),TRGTBL(1501)
      DIMENSION FMTMAP(7),FMTIND(5)
      INTEGER FMTMAP,FMTIND,BETAX,BETAZ
      DIMENSION RR(101),IX2(85)
      DIMENSION IX1(85),WORDS(24)
      DIMENSION SS(1996),NSS(1999)
      DIMENSION ALPHA(2,12)
      DIMENSION NR(85,47)
      COMMON TRGTBL,NXIS,NXX,NYX,NZ,XMIN,XMAX,YMIN,YMAX,ZMIN
      COMMON ZMAX,DELTAX,DY,DZ,AMAX,BMAX,JMAX,KMAX,KMIN,NMEMMX
      COMMON NEGINA,NCD,NEF,NEFTBL,NX,NY,NCUE,N1,RMAX,NLAST,NFPLT
      COMMON YFL,N11,NN1,III,IVXX,IXX,NN4,NN5,NN6,ISST,AXXX,BXXX
      COMMON CXXX,COSA,COSB,COSC,NSEC,NGRID,NSKIP,CNTRNO,CNTNEG
      COMMON NEXPAN,WORDS,NR,X,Y,RR
      character*28 FMTSFY
      equivalence ( FMTSFY,FMTMAP )
      character    LINSFY(130)
      character*30 FMT001
      character    ALPHASFY(2,12)
      EQUIVALENCE (NR(1),SS(1)),(NR(1997),NSS(1))
      INTEGER ALPHA1,ALPHA2,ALPHA3,BETAX2,BETAX1,
     1 BETAX3,BETAZ1,BETAZ2,BETAZ3
      INTEGER ALPHAX,ALPHAY,ALPHAZ,CHAR,ALPHA
Csfy  DATA FMTMAP/'(1H ',',00X',',2HY','F8.5'
Csfy 1 ,' ,3X',',00A','1)  '/
      DATA FMTMAP/'(1H ',',00X',',1HY',',F8.'
     1 ,'5,3X',',00A','1)  '/
      DATA FMTIND/'(1H ',',00X',',13X','00A '
     1 ,'1)  '/
      DATA ALPHAX/' X='/
      DATA ALPHAY/' Y='/
      DATA ALPHAZ/' Z='/
      DATA ALPHA/'1','A','2','B','3','C','4','D','5','E','6','F','7','G'
     1,'8','H','9','I','0','O','+','-',' ',' '/
      DATA ALPHASFY / '1','A',
     &                '2','B',
     &                '3','C',
     &                '4','D',
     &                '5','E',
     &                '6','F',
     &                '7','G',
     &                '8','H',
     &                '9','I',
     &                '0','O',
     &                '+','-',
     &                ' ',' ' /
      DATA BETAX/1HX/
      DATA BETAZ/1HZ/
      DATA ALPHA1/2H 1,2H A,2H 2,2H B,2H 3,2H C,2H 4,2H D,2H 5,2H E,2H 6
     2,2H F,2H 7,2H G,2H 8,2H H,2H 9,2H I,2H 0,2H O,2H +,2H -,2H  ,2H  /
      DATA ALPHA2/3H  1,3H  A,3H  2,3H  B,3H  3,3H  C,3H  4,3H  D,3H  5,
     33H  E,3H  6,3H  F,3H  7,3H  G,3H  8,3H  H,3H  9,3H  I,3H  0,3H  O,
     43H  +,3H  -,3H   ,3H   /
      DATA ALPHA3/4H   1,4H   A,4H   2,4H   B,4H   3,4H   C,4H   4,4H
     5D,4H   5,4H   E,4H   6,4H   F,4H   7,4H   G,4H   8,4H   H,4H   9,
     64H   I,4H   0,4H   O,4H   +,4H   -,4H    ,4H    /
      DATA BETAX1/2H X/
      DATA BETAX2/3H  X/
      DATA BETAX3/4H   X/
      DATA BETAZ1/2H Z/
      DATA BETAZ2/3H  Z/
      DATA BETAZ3/4H   Z/
C
3001  FORMAT (69H                                                     ++
     1++++++++++++++)
  300 FORMAT (6H RMAX=E13.6,35H                                  '2XA4,F
     18.5,36H'                                C =E13.6)
2125  FORMAT (1H )
 2225 FORMAT (1H1)
 3012 FORMAT (14H              29(1XI3))
C
C
      READ(60,800)NHORI,NVERT,NDISP,CMIN,CINT
  800 FORMAT(///3I5,2F5.0)
      IF(NHORI.GT.3) NHORI=3
      WRITE(61,900)
  900 FORMAT(//1H 'ALPHAMERIC FOURIER PLOT'//1H 'POSITIVE POINTS NUMERIC
     1    0 + 1 2 3 4 5 6 7 8 9'//1H 'NEGATIVE POINTS ALPHABETIC O - A B
     1 C D E F G H I'//)
      WRITE(61,901)NHORI,NVERT,NDISP,CMIN,CINT
  901 FORMAT(' NHORI=',I3,'   NVERT=',I3,'   NDISP=',I3,'   CMIN=',
     1 F6.0,'   CINT=',F6.0/)
C
C  SET UP ARRAYS FOR A1,A2,A3 OR A4
C
      IF(NHORI.EQ.0) GO TO 21337
      IF(NHORI -2) 721,722,723
  721 DO 1721 I=1,12
      DO 1721 J=1,2
 1721 ALPHA(J,I)=ALPHA1(J,I)
      BETAX=BETAX1
      BETAZ=BETAZ1
      GO TO 21337
  722 DO 1722 I=1,12
      DO 1722 J=1,2
 1722 ALPHA(J,I)=ALPHA2(J,I)
      BETAX=BETAX2
      BETAZ=BETAZ2
      GO TO 21337
  723 DO 1723 I=1,12
      DO 1723 J=1,2
 1723 ALPHA(J,I)=ALPHA3(J,I)
      BETAX=BETAX3
      BETAZ=BETAZ3
C
C
21337 NPRT=10*ISST(9)+ISST(10)
      IF (NPRT)4,4,5
    4 NPRT=118/(1+NHORI)
      IF(NHORI.EQ.0) NPRT=85
      GO TO 5510
    5 IF(NPRT.GT.(118/(1+NHORI)))GO TO 4
 5510 IF (AMAX+BMAX)5999,1901,5999
 1901 AMAX=1000.
 5999 REWIND IVXX
  196 IF (RMAX)202,200,202
  200 READ(IVXX         ) CODE,RMAX
     1     ,     NX,NY,XMIN,XMAX,DELTAX,ZMAX,Z,NXIS,YMIN,YMAX,DY,AMAX,
     1BMAX,DZ,ZMIN,NZ
      IF(CODE.LT.0.0) GO TO 202
      DO 201 L=1,NY
  201 READ(IVXX)
      GO TO 200
  202 C=999./RMAX
  310 REWIND IVXX
      IF(NXIS.NE.10) GO TO 22834
Csfy  FMTMAP(3)=FMTMAP(3)+16777216
C22834 FMTMAP(7)=FMTMAP(7)+NHORI*4398046511104D00
Csfy 22834 FMTMAP(7)=FMTMAP(7)+NHORI*1
22834 continue
      FMTIND(5)=FMTMAP(7)
      NTEMP0=NPRT/10
C     FMTMAP(6)=FMTMAP(6)+NTEMP0*68719476736D00
      FMTMAP(6)=FMTMAP(6)+NTEMP0*256+(NPRT-NTEMP0*10)*65536
C    A+(NPRT-NTEMP0*10)*1073741824
      FMTIND(4)=FMTMAP(6)
22835 NBLANK=1
      IF(NDISP.LE.0) NBLANK=118-(NHORI+1)*NPRT
      IF(((NHORI+1)*NPRT+IABS(NDISP)*NY).GT.118) NDISP=0
      NTEMP3=NBLANK/10
      NTEMP4=NBLANK-NTEMP3*10
C     FMTIND(2)=FMTIND(2)+NTEMP3*68719476736D00+NTEMP4*1073741824
      FMTIND(2)=FMTIND(2)+NTEMP3*256 + NTEMP4*65536
C
C  FINISH OF INITIALISATION
C
 2281 READ(IVXX) CODE,CODD,
     1           NX,NY,XMIN,XMAX,DELTAX,ZMAX,Z,NXIS,YMIN,YMAX,DY,AMAX,
     1BMAX,DZ,ZMIN,NZ
      IF(CODE.LT.0.0) GO TO 219
  404 NVAR=0
  405 X(1)=XMIN
      DO 204 J=2,NX
204   X(J)=X(J-1)+DELTAX
      WRITE (III,2225)
      WRITE (III,3001)
      IF (NXIS-10)2041,2042,2043
 2041 CHAR=ALPHAZ
      GO TO 2044
 2042 CHAR=ALPHAY
      GO TO 2044
 2043 CHAR=ALPHAX
 2044 WRITE (III,300)RMAX,CHAR,Z,C
      WRITE (III,3001)
Csfy  PRINT 2125
 207  NP=1
      NQ=NPRT
 216  IF (NX-NQ)98,218,218
98    NQ=NX
C
C  REINITIALISE FORMATS
C
  218 L=NQ-NP+1
      NTEMP0=NTEMP3
      NTEMP1=NTEMP4
      FMTMAP(2)=FMTIND(2)
      FMTIND(4)=FMTMAP(6)
C
C
      IF(NXIS.GT.10) GO TO 2182
      CHAR=BETAX
      GO TO 2183
 2182 CHAR=BETAZ
 2183 DO 22182 I=1,L
22182 IX1(I)=CHAR
      WRITE(III,FMTIND) (IX1(I),I=1,L)
      K=0
      DO 3 J=NP,NQ
      K=K+1
      IF(X(J).LT.0.0) GO TO 333
      IX1(K)=ALPHA(1,11)
      GO TO 3
  333 IX1(K)=ALPHA(2,11)
    3 CONTINUE
      WRITE(III,FMTIND) (IX1(K),K=1,L)
C     FMTIND(4)=FMTIND(4)+134217728
      FMTIND(4)=FMTIND(4)+134217728
      K=0
      DO 1333 I=NP,NQ
      K=K+1
      IX2(K)=     AINT(X(I)*100.0+SIGN(0.1,X(I)))
 1333 IX2(K)=IABS(IX2(K))
      DO 2333 I=1,L
      IX1(I)=IX2(I)/10
 2333 IX2(I)=IX2(I)-10*IX1(I)
      WRITE(III,FMTIND) (IX1(I),I=1,L)
      WRITE(III,FMTIND) (IX2(I),I=1,L)
      WRITE(III,2125)
C
C     READ ANSWERS 1 LINE AT A TIME
C
  213 DO 2123 N=1,NY
  211 READ(IVXX) CODE,CODE,Y(1),(SS(J),J=1,NX)
C
C     FIND ALPHAMERIC SYMBOL FOR EACH ANSWER ON THE LINE
C
      DO 212 J=NP,NQ
      SS(J)=C*SS(J)
      IND=1
      IF(SS(J).LT.CMIN)GO TO 735
      IF(SS(J).GE.0.)GO TO 710
      IND=2
      SS(J)=-SS(J)
  710 DO 720 L=1,9
      IF(SS(J).GE.((L*100)+CINT-1) .AND.SS(J).LT.((L+1)*100))  GO TO 740
  720 CONTINUE
      IF(SS(J).GE.(CINT-1).AND.SS(J).LT.100) GO TO 730
      IF(SS(J).LE.(50.-CINT/2.))GO TO 725
      GO TO 735
  730 L=11
      GO TO 740
  725 L=10
      GO TO 740
  735 L=12
  740 NSS(J)=ALPHA(IND,L)
      JSFY = J - NP + 1
      LINSFY(JSFY) = ALPHASFY(IND,L)
  212 CONTINUE
C
C     OUTPUT 1 LINE OF ALPHAMERIC PLOT ON THE LINE PRINTER
C
      write ( FMT001,9001 ) NBLANK,NPRT
 9001 format ( '(1x,',i2,'x,1hY,f8.5,3x,',i2,'a1)' )
Csfy  WRITE(III,FMTMAP) Y(1),(NSS(J),J=NP,NQ)
      WRITE(III,FMT001) Y(1),(LINSFY(J),J=1,NPRT)
      IF(NVERT.EQ.0) GO TO 22936
      DO 22836 J=1,NVERT
Csfy  PRINT 2125
22836 CONTINUE
22936 IF(NDISP.EQ.0) GO TO 2123
      NTEMP2=NTEMP1+NDISP
      IF(NTEMP2.LT.0) GO TO 2111
      IF(NTEMP2.LE.9) GO TO 8001
      NTEMP2=NTEMP2-10
 2111 NTEMP0=NTEMP0+ISIGN(1,NDISP)
C     FMTMAP(2)=FMTMAP(2)+68719476736*ISIGN(1,NDISP)
      FMTMAP(2)=FMTMAP(2)+256*ISIGN(1,NDISP)
      IF(NTEMP2.LT.0) NTEMP2=10+NTEMP2
C     FMTMAP(2)=FMTMAP(2)+1073741824*(NTEMP2-NTEMP1)
      FMTMAP(2)=FMTMAP(2)+65536*(NTEMP2-NTEMP1)
      NTEMP1=NTEMP2
      GO TO 2123
 8001 NTEMP1=NTEMP1+NDISP
C     FMTMAP(2)=FMTMAP(2)+1073741824*NDISP
      FMTMAP(2)=FMTMAP(2)+65536*NDISP
 2123 CONTINUE
C
C     LOOP BACK FOR ANOTHER LINE
C
C
C     MORE THAN  29 ANSWERS/LINE PLOT THE REMAINDER
C
 214  IF (NX-NQ)2281,2281,217
  217 WRITE (III,2225)
      NBSP=NY
      DO 220 M=1,NBSP
  220 BACKSPACE IVXX
      NP=NP+NPRT
      NQ=NQ+NPRT
      GO TO 216
  219 REWIND IVXX
      RETURN
      END


      SUBROUTINE IPEAK(U)
      COMMON /PEAK/ RMIN,NPEAK,DX,DY,DZ
      COMMON DUMMY(1533),IN,III,IVXX,IXX,DUM(3),ISST(10),DM(12),WORDS(24
     A ),IBOR,NNZ,RSAVE
      DIMENSION LABEL(100),LNO(100),XX(100),YY(100),ZZ(100),R(100),
     A RSAVE(2700)     ,XXINT(3,100),IBOR(100),BORDER(3),U(3,3),XO(3)
      INTEGER NNX,NNY,NNZ
      DIMENSION NNX(100),NNY(100),NNZ(100)
      EQUIVALENCE (DUMMY(1),LABEL(1)),(DUMMY(101),LNO(1)),(DUMMY(201),
     A XX(1)),(DUMMY(301),YY(1)),(DUMMY(401),ZZ(1)),(DUMMY(501),R(1)),
     B (DUMMY(601),XXINT(1)),(DUMMY(901),NNX(1)),(DUMMY(1001),NNY(1))
       DATA BORDER/1H ,2H +,4H +++/
C
      IF(NPEAK.EQ.0) GO TO 880
      IF(NPEAK.GT.100) NPEAK=100
      REWIND IVXX
      REWIND IXX
      DO 1 I=1,NPEAK
    1 READ(IXX) LABEL(I),LNO(I),XX(I),YY(I),ZZ(I),IBOR(I),
     A                                            NNX(I),NNY(I),NNZ(I)
      NZ = 0
   14 READ(IVXX) CODE,CODD,NX,NY,CODD,CODD,DX,CODD,CODD,NXIS,CODD,
     A  CODD,DY,CODD,CODD,DZ
      NZ=NZ+1
      NUM=0
      IF(CODE.LT.0.0) GO TO 700
  100 NUM=NUM+1
      IF(NUM.GT.NY) GO TO 14
      READ(IVXX) CODE,CODE,Y,(R(I),I=1,NX)
      DO 21 I=1,NPEAK
      IF(IBOR(I).EQ.2) GO TO 21
      J=NNZ(I)
      IF(IABS(J-NZ).GT.1) GO TO 21
      J=NZ-J+2
      K=NNY(I)-1
      IF(IABS(K-NUM).GT.1) GO TO 21
      K=NUM-K+2
      M=NNX(I)-2
      DO 20 L=1,3
      KKK=-39+27*I+9*J+3*K+L
20    RSAVE(KKK)=R(M+L)
   21 CONTINUE
      GO TO 100
C
  700 REWIND IVXX
      REWIND IXX
      IF(NXIS.EQ.1) GO TO 200
      DO 201 I=1,NPEAK
      IF(NXIS.EQ.10) GO TO 202
      CALL ROTATE(XX(I),ZZ(I))
      GO TO 201
  202 CALL ROTATE(YY(I),ZZ(I))
  201 CONTINUE
  200 WRITE(III,204)
  210 WRITE(III,211)
  240 IF(ISST(2).EQ.0) GO TO 1240
      WRITE(III,2240)
      GO TO 3250
 1240 WRITE(III,1260)
 3250 DO1250 I=1,NPEAK
      IF(IBOR(I).EQ.1) GO TO 2250
  400 XXINT(1,I)=XX(I)
      XXINT(2,I)=YY(I)
      XXINT(3,I)=ZZ(I)
      GO TO 250
 2250 CALL EDINT(XX(I),YY(I),ZZ(I),RSAVE(27*I-26),XXINT(1,I),IBOR(I))
      IF(IBOR(I).EQ.3) GO TO 400
  250 IF(NXIS-10) 470,480,490
  490 CALL ROTATE(XXINT(1,I),XXINT(3,I))
      GO TO 470
  480 CALL ROTATE(XXINT(2,I),XXINT(3,I))
      IDORI=IBOR(I)
  470 WRITE(III,251)I,(XXINT(J,I),J=1,3),BORDER(IDORI)
      DO 510 J=1,3
      XO(J)=0.0
      DO 510 K=1,3
  510 XO(J)=XO(J)+U(K,J)*XXINT(K,I)
      IF(ISST(2).EQ.0) GO TO 610
      DPAT=0.0
      GO TO 1250
  610 DPAT=SQRT(XO(1)*XO(1)+XO(2)*XO(2)+XO(3)*XO(3))
 1250 WRITE(III,1270) XO,DPAT
      NPCH=4*ISST(7)
C
  204 FORMAT(///'1INTERPOLATED COORDINATES FOR PEAKS (THOSE MARKED WITH
     AASTERISK(S) HAVE NOT BEEN INTERPOLATED)'//)
  211 FORMAT('0PEAK NO.'13X'X'12X'Y'12X'Z')
  251 FORMAT('0'I5,9X3(2XF9.4,2X),A4)
  252 FORMAT(2A4,14X,'1.000000'3F10.6)
 1260 FORMAT('+'77X'X(A)'4X'Y(A)'4X'Z(A)'8X'D(000)')
 2240 FORMAT('+'77X'OLDX'4X'OLDY'4X'OLDZ')
 1270 FORMAT('+'73X3F8.3,5XF8.3)
C
  300 DO 320 I=1,NPEAK
      IF(NPCH.LE.0) GO TO 320
      NPCH=NPCH-1
C     PUNCH  252,LABEL(I),LNO(I),(XXINT(J,I),J=1,3)
  320 WRITE(IXX) LABEL(I),LNO(I),(XXINT(J,I),J=1,3)
      REWIND IXX
  880 RETURN
      END


      SUBROUTINE FRPRNT
      DIMENSION X( 103),ISST(10)
      DIMENSION Y(101),TRGTBL(1501)
      DIMENSION RR(101),IX2(29)
      DIMENSION IX1(29),WORDS(24)
      DIMENSION SS(1996),NSS(1999)
      DIMENSION NR(85,47)
      COMMON TRGTBL,NXIS,NXX,NYX,NZ,XMIN,XMAX,YMIN,YMAX,ZMIN
      COMMON ZMAX,DELTAX,DY,DZ,AMAX,BMAX,JMAX,KMAX,KMIN,NMEMMX
      COMMON NEGINA,NCD,NEF,NEFTBL,NX,NY,NCUE,N1,RMAX,NLAST,NFPLT
      COMMON YFL,N11,NN1,III,IVXX,IXX,NN4,NN5,NN6,ISST,AXXX,BXXX
      COMMON CXXX,COSA,COSB,COSC,NSEC,NGRID,NSKIP,CNTRNO,CNTNEG
      COMMON NEXPAN,WORDS,NR,X,Y,RR
      EQUIVALENCE (NR(1),SS(1)),(NR(1997),NSS(1))
      INTEGER ALPHAX,ALPHAY,ALPHAZ,CHAR
      DATA ALPHAX/' X='/
      DATA ALPHAY/' Y='/
      DATA ALPHAZ/' Z='/
  303 FORMAT (14H              29A4)
  302 FORMAT (3H Y=F8.5,3H   29I4)
 3021 FORMAT (3H Z=F8.5,3H   29I4)
3001  FORMAT (69H                                                    ++
     1+++++++++++++*)
  300 FORMAT (6H RMAX=E13.6,35H                                  +2XA4,F
     18.5,36H+                                C =E13.6)
2125  FORMAT (1H )
 2225 FORMAT (1H1)
 3012 FORMAT (14H              29(1XI3))
      NPRT=10*ISST(9)+ISST(10)
      IF (NPRT)4,4,5
 4    NPRT=29
      GO TO 6
 5    IF (NPRT-29)6,6,4
 6    IF (ISST(6)-2)1,2,1
 1    ISZ=1
      GO TO 5510
 2    ISZ=2
 5510 IF (AMAX+BMAX)5999,1901,5999
 1901 AMAX=1000.
 5999 REWIND IVXX
  196 IF (RMAX)202,200,202
  200 READ (IVXX)CODE,RMAX
     1     ,     NX,NY,XMIN,XMAX,DELTAX,ZMAX,Z,NXIS,YMIN,YMAX,DY,AMAX,
     1BMAX,DZ,ZMIN,NZ
      IF(CODE.LT.0.0) GO TO 202
      DO 201 L=1,NY
  201 READ(IVXX)
      GO TO 200
  202 C=999./RMAX
  310 REWIND IVXX
 2281 READ(IVXX) CODE,CODD,
     1           NX,NY,XMIN,XMAX,DELTAX,ZMAX,Z,NXIS,YMIN,YMAX,DY,AMAX,
     1BMAX,DZ,ZMIN,NZ
      IF(CODE.LT.0.0) GO TO 219
      IF (NX- 85)6023,6023,404
 6023 IF (NY- 47)400,400,404
  404 NVAR=0
      GO TO 405
  400 NVAR=1
      DO 440 I=1,NY
      READ(IVXX) CODE,CODE,Y(I),(RR(J),J=1,NX)
      DO 440 J=1,NX
      IF (RR(J))438,439,439
 438  NR(J,I)=C*RR(J)-0.5
      GO TO 440
 439  NR(J,I)=C*RR(J)+0.5
  440 CONTINUE
  405 X(1)=XMIN
      DO 204 J=2,NX
204   X(J)=X(J-1)+DELTAX
      WRITE (III,2225)
      WRITE (III,3001)
      IF (NXIS-10)2041,2042,2043
 2041 CHAR=ALPHAZ
      GO TO 2044
 2042 CHAR=ALPHAY
      GO TO 2044
 2043 CHAR=ALPHAX
 2044 WRITE (III,300)RMAX,CHAR,Z,C
      WRITE (III,3001)
 207  NP=1
      NQ=NPRT
 216  IF (NX-NQ)98,218,218
98    NQ=NX
218   IF (NXIS-10)2181,2181,2182
 2181 CHAR=ALPHAX
      GO TO 2183
 2182 CHAR=ALPHAZ
 2183 WRITE (III,303)CHAR,CHAR,CHAR,CHAR,CHAR,CHAR,CHAR,CHAR,CHAR,CHAR,
     1CHAR,CHAR,CHAR,CHAR,CHAR,CHAR,CHAR,CHAR,CHAR,CHAR,CHAR,CHAR,CHAR,
     1CHAR,CHAR,CHAR,CHAR,CHAR,CHAR
      DO 3 J=NP,NQ
      K=1+MOD(J-1,NPRT)
      IX1(K)=AINT(100.0*X(J)+SIGN(0.1,X(J)))
 3    IX2(K)=IFIX(100000.*ABS(X(J))+.1)-1000*IABS(IX1(K))
      NNP=1+MOD(NP-1,NPRT)
      NNQ=1+MOD(NQ-1,NPRT)
      WRITE (III,3012)(IX1(K),K=NNP,NNQ)
      WRITE (III,3012)(IX2(K),K=NNP,NNQ)
      WRITE (III,2125)
      IF (NVAR)213,213,450
  450 IF (NXIS-10)459,469,459
 459  IF (ISZ-1)460,460,486
 486  DO 487 I=1,NY
      WRITE (III,302)Y(I),(NR(J,I),J=NP,NQ)
 487  WRITE (III,2125)
      GO TO 480
  460 DO 462 I=1,NY
  462 WRITE (III,302)Y(I),(NR(J,I),J=NP,NQ)
      GO TO 480
 469  IF (ISZ-1)470,470,488
 488  DO 489 I=1,NY
      WRITE (III,3021)Y(I),(NR(J,I),J=NP,NQ)
 489  WRITE (III,2125)
      GO TO 480
  470 DO 472 I=1,NY
  472 WRITE (III,3021)Y(I),(NR(J,I),J=NP,NQ)
 480  IF (NX-NQ)2281,2281,485
  485 WRITE (III,2225)
      NP=NP+NPRT
      NQ=NQ+NPRT
      GO TO 216
  213 DO 2123 N=1,NY
  211 READ(IVXX) CODE,CODE,Y(1),(SS(J),J=1,NX)
      DO 212 J=NP,NQ
      IF (SS(J))230,231,231
 230  NSS(J)=C*SS(J)-0.5
      GO TO 212
 231  NSS(J)=C*SS(J)+0.5
  212 CONTINUE
      IF (NXIS-10)2111,2112,2111
 2111 IF (ISZ-1)1000,1000,1001
 1000 WRITE (III,302)Y(1),(NSS(J),J=NP,NQ)
      GO TO 2123
 1001 WRITE (III,302)Y(1),(NSS(J),J=NP,NQ)
      WRITE (III,2125)
      GO TO 2123
 2112 IF (ISZ-1)1002,1002,1003
 1002 WRITE (III,3021)Y(1),(NSS(J),J=NP,NQ)
      GO TO 2123
 1003 WRITE (III,3021)Y(1),(NSS(J),J=NP,NQ)
      WRITE (III,2125)
 2123 CONTINUE
 214  IF (NX-NQ)2281,2281,217
  217 WRITE (III,2225)
      NBSP=NY
      DO 220 M=1,NBSP
  220 BACKSPACE IVXX
      NP=NP+NPRT
      NQ=NQ+NPRT
      GO TO 216
  219 REWIND IVXX
      RETURN
      END


      SUBROUTINE GEMI(A,B,X,N,NEROR)
      REAL A,B
      DIMENSION A(20,20),B(3)
      REAL X(3)
      DO 10 I=1,3
   10 A(I,4)=B(I)
      DO 1 III=1,3
      A(4,4)=1.D0/A(1,1)
      DO 2 I=1,3
    2 A(4,I)=A(1,I+1)*A(4,4)
      DO 3 I=1,2
      DO 4 J=1,4
    4 A(I,J)=-A(I+1,1)*A(4,J)
      DO 3 J=1,3
    3 A(I,J)=A(I,J)+A(I+1,J+1)
      DO 1 J=1,4
    1 A(3,J)=A(4,J)
      DO 5 I=1,3
    5 X(I)=A(I,1)
      RETURN
      END


      SUBROUTINE GP1(PLANE,D,AA,U)
Csfy  DIMENSION PLANE(3),A(6)
      DIMENSION PLANE(3)
      DIMENSION U1(3),U2(3)
      DIMENSION U3(3),U(3,3)
      DIMENSION AA(3,3),BB(3,3)
C     SUBROUTINE TO FIND TRANSFORMATION
C     MATRIX TO AXES APPROPRIATE TO PLANE
C     CELL2=U*CELL1
C     CELL2(1)=NORMAL TO PLANE
C     SET UP METRICS
      CALL MI3(AA,BB,DETA)
      IF (D)4,3,3
    4 DO 5 I=1,3
    5 PLANE(I)=-PLANE(I)
      D=-D
    3 ALPHA=1./SQRT(VMV(PLANE,BB,PLANE))
      DO 6 I=1,3
    6 PLANE(I)=ALPHA*PLANE(I)
      D=ALPHA*D
C     PLANE IS NOW NORMALIZED WITH POSITIVE D
      CALL VM(PLANE,BB,U1)
C     U1 IS UNIT VECTOR PERPENDICULAR TO PLANE
      K=1
      IF (ABS(PLANE(1))-ABS(PLANE(2)))7,8,8
    8 K=2
    7 IF (ABS(PLANE(K))-ABS(PLANE(3)))9,10,10
   10 K=3
    9 FACTOR=-SIGN(PLANE(K),U1(K))
      DO 11 I=1,3
      IF (I-K)12,13,12
   13 U2(I)=1.+FACTOR*U1(I)
      GO TO 11
   12 U2(I)=FACTOR*U1(I)
   11 CONTINUE
      CALL NORM(AA,U2)
C     WE NOW HAVE NORMALIZED SECOND VECTOR
      CALL CROSS(BB,U1,U2,U3)
      CALL NORM(AA,U3)
C     WE NOW HAVE AN ORTHONORMAL SET
      DO 1 I=1,3
      U(1,I)=U1(I)
      U(2,I)=U2(I)
    1 U(3,I)=U3(I)
      RETURN
      END


      SUBROUTINE MI3(A,AI,D)
      DIMENSION A(3,3),AI(3,3)
      DO 3 I=1,3
      K=1+MOD(I+1,3)
      M=1+MOD(I+3,3)
      DO 3 J=1,3
      L=1+MOD(J+1,3)
      N=1+MOD(J+3,3)
    3 AI(I,J)=A(L,K)*A(N,M)-A(L,M)*A(N,K)
      D=0.0
      DO 1 I=1,3
    1 D=D+A(I,1)*AI(1,I)
      DO 2 I=1,3
      DO 2 J=1,3
    2 AI(I,J)=AI(I,J)/D
      RETURN
      END


      SUBROUTINE MV(A,B,C)
      DIMENSION B(3),C(3)
      DIMENSION A(3,3)
      DO 2 I=1,3
      C(I)=0
      DO 2 J=1,3
    2 C(I)=A(I,J)*B(J)+C(I)
      RETURN
      END


      SUBROUTINE NORM(AA,U)
      DIMENSION U(3),AA(3,3)
      ALPHA=SQRT(1./VMV(U,AA,U))
      DO 1 I=1,3
    1 U(I)=ALPHA*U(I)
      RETURN
      END


      SUBROUTINE RIP(U,VOL)
      DIMENSION MTX(24),MTY(24)
      DIMENSION MTZ(24),J(24)
      DIMENSION K(24),L(24)
      DIMENSION A(24),B(24)
      DIMENSION RN(3),WORDS(24)
      DIMENSION JX(24),KX(24)
      DIMENSION LX(24),AX(24)
      DIMENSION BX(24),TRGTBL(1501)
      DIMENSION ISST(10),COUGH(9)
      DIMENSION M1(24,3),M2(24,3)
      DIMENSION M3(24,3),U(3,3)
      DIMENSION G(3,3),WORD1(4),WORD2(4)
      DIMENSION HJH(10),HJK(10),HJL(10),SINTOL(10),FO(10),SFCA(10),
     1 COSPHA(10),SINPHA(10),DELTA(10),JUNK(10),JANK(10)
      INTEGER   HJH,HJK,HJL,DELTA,JANK
      INTEGER JHALF,KHALF,LHALF
      DIMENSION JHALF(24),KHALF(24),LHALF(24)
      COMMON TRGTBL,NORIEN,NXX,NYX,NZ,XMIN,XMAX,YMIN,YMAX,ZMIN
      COMMON ZMAX,DX,DY,DZ,AMAX,BMAX,JMAX,KMAX,KMIN,NMEMMX,NEGINA
      COMMON NCD,NEF,NEFTBL,NX,NY,NCUE,N1,RMAX,NLAST,NFPLT,YFL
      COMMON N11,IN,III,IVXX,IXX,NN5,NN6,NDATAP,ISST,AXX,BXX
      COMMON C,COSA,COSB,COSC,NSEC,NGRID,NSKIP,CNTRNO,CNTNEG
      COMMON NEXPAN,WORDS,MTX,MTY,MTZ,J,K,L,A,B,JX,KX,LX,AX,BX
      COMMON M1,M2,M3,RN,RD,G,COUGH
      COMMON FO,JUNK,SFCA,SINTOL,COSPHA,SINPHA,HJH,HJK,HJL,DELTA,JANK
      COMMON JHALF,KHALF,LHALF     ,AHALF(24),BHALF(24)
      DATA WORD1/'MAGN','ETIC','  FI','LE'/
      DATA WORD2/'CARD','S','  ','  '/
 80   FORMAT (5I5,F5.2)
   88 FORMAT (20H0INPUT DATA IS FROM 4A4)
   90 FORMAT (9F5.5,5X,I3)
   91 FORMAT (13H0X GOES FROM F8.5,4H TO F8.5,13H IN STEPS OF F8.5    /
     113H0Y GOES FROM F8.5,4H TO F8.5,13H IN STEPS OF F8.5    /13H0Z GOE
     1S FROM F8.5,4H TO F8.5,13H IN STEPS OF F8.5/
     116H0ORIENTATION IS I3)
 109  FORMAT (3(F15.6,3I3))
 3008 FORMAT (50X,45H(THE INDICES ARE TRANSFORMED BY THE U MATRIX))
  115 FORMAT (39H0     H     K     L      A         B   /)
  116 FORMAT (39H0     H     L     K      A         B   /)
  117 FORMAT (39H0     L     K     H      A         B   /)
  203 FORMAT(3I5,5X,2F10.5)
  305 FORMAT (1H ,3I6,2F10.2)
  306 FORMAT (1H )
 2000 FORMAT (31H0NO INPUT DATA = NO CALCULATION)
 2001 FORMAT (43H0THERE ARE NO EXPANDED DATA END CALCULATION)
  352 FORMAT (5H0THE I5,47H INDEPENDENT REFLECTIONS HAVE BEEN EXPANDED T
     1O I6,8H TERMS. )
  355 FORMAT (23H0THE CD TABLE REQUIRES I5,14H MEMORY UNITS.)
  357 FORMAT (23H0THE EF TABLE REQUIRES I5,14H MEMORY UNITS.)
  361 FORMAT (36H0THE MAXIMUM SIZE OF THE EF TABLE ISI5,57H MEMORY UNITS
     1.  THE PROGRAM WILL SPLIT EACH SECTION INTO I2,7H PARTS.)
 4008 FORMAT (47H0INSUFFICIENT STORAGE AVAILABLE. TERMINATE JOB.)
      READ (IN,80)NEQUIV,NREF,NTAPIN,NONTAP,NREJ,SNLMAX
      IF(SNLMAX.EQ.0.) SNLMAX=3.0
      WRITE(III,980)NEQUIV,NREF,NTAPIN,NONTAP,NREJ,SNLMAX
  980 FORMAT('0NEQUIV=',I3,3X'NREF=',I4,3X'NTAPIN=',I4,3X'NONTAP=',
     A I5,3X'NREJ=',I3,3X'SNLMAX=',F6.2)
      NTAPSV=NTAPIN
      IF (NTAPIN-100)3002,3002,3000
 3000 NTAPIN=NTAPIN-100
      GO TO 3002
 3002 IF (NTAPIN-10)18,16,16
   16 IF (NONTAP)18,17,18
   17 NONTAP=10000
   18 REWIND IVXX
      NEXPAN=0
      NTOTAL=0
      NDATAP=2
      NUMB=10
      RESID=0.0
      NHOLD=0
      ASSIGN 311 TO NXIT
      NNUMB=0
      AMAX=0.
      BMAX=0.
      NZERO=0
      JMAX=0
      KMAX=0
      KMIN=0
   81 IF (NTAPIN-1)85,84,84
84    WRITE(III,88)WORD1
      GO TO 87
85    WRITE(III,88)WORD2
87    IF (ISST(2)-1)3003,3003,3004
 3003 READ (IN,90)XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,DX,DY,DZ,NORIEN
      GO TO 3005
 3004 READ (IN,90)(COUGH(I),I=1,9),IHACK
      NORIEN=100
 3005 CALL CHEKUP(XMIN,XMAX,NORIEN,100,DX)
      CALL CHEKUP(YMIN,YMAX,NORIEN,10,DY)
      CALL CHEKUP(ZMIN,ZMAX,NORIEN,1,DZ)
  104 WRITE (III,91)XMIN,XMAX,DX,YMIN,YMAX,DY,ZMIN,ZMAX,DZ,NORIEN
      IF (NORIEN-10)98,96,94
   94 CALL ROTATE(XMIN,ZMIN)
      CALL ROTATE(XMAX,ZMAX)
      CALL ROTATE(DX,DZ)
      GO TO 98
   96 CALL ROTATE(YMIN,ZMIN)
      CALL ROTATE(YMAX,ZMAX)
      CALL ROTATE(DY,DZ)
   98 NX=(XMAX-XMIN)/DX+1.001
      NY=(YMAX-YMIN)/DY+1.001
      NZ=(ZMAX-ZMIN)/DZ+1.001
  122 WRITE(III,27918)
27918 FORMAT('0SYMMETRY INFORMATION'/)
      DO 139 I=1,NEQUIV
      READ (IN,109,END=10)  (RN(N),M1(I,N),M2(I,N),M3(I,N),N=1,3)
C     IF(EOFCKF(IN).EQ.1)STOP
6128  WRITE (III,109)(RN(N),M1(I,N),M2(I,N),M3(I,N),N=1,3)
      IF (M1(I,1)+M2(I,1))170,400,170
 400  IF (M1(I,1))171,179,171
  170 IF (M1(I,2)+M2(I,2))179,401,179
 401  IF (M1(I,2))171,179,171
C     HEXAGONAL SITUATIONS
  171 IF (M1(I,1)+M2(I,2))173,172,173
  172 NSAVE=M2(I,1)
      M2(I,1)=M1(I,2)
      M1(I,2)=NSAVE
      GO TO 179
  173 NSAVE=M2(I,2)
      M2(I,2)=M1(I,1)
      M1(I,1)=NSAVE
  179 DO 129JJ=1,3
  132 IF (RN(JJ))133,129,129
  133 RN(JJ)=1.+RN(JJ)
      GO TO 132
  129 CONTINUE
      MTX(I)=12.*RN(1)+.1
      MTY(I)=12.*RN(2)+.1
  139 MTZ(I)=12.*RN(3)+.1
  110 IF (NREF)111,201,111
  111 IF (NORIEN-10)112,113,114
  112 WRITE (III,115)
      GO TO 3006
  113 WRITE (III,116)
      GO TO 3006
  114 WRITE (III,117)
 3006 IF (ISST(2)-1)201,3007,3007
 3007 WRITE (III,3008)
  201 IF(NTAPIN-10)202,180,180
  180 IF (NONTAP)320,320,181
  181 NONTAP=NONTAP-1
C     IF(NUMB.LT.10) GO TO 86425
C86426 NUMB=0
C
C     THIS STATEMENT WILL READ BOTH AN IED FILE AND A DP3(OR EQUIVALENT)
C
C     READ(NDATAP)           (FO(MN),JUNK(MN),SFCA(MN),SINTOL(MN),
C    1 COSPHA(MN),SINPHA(MN),MN=1,10), (HJH(MON),HJK(MON),HJL(MON),
C    2 DELTA(MON),JANK(MON),MON=1,10)
C     IF(EOFCKF(NDATAP).EQ.1)GOTO33986
C86425 NUMB=NUMB+1
C     JJ=HJH(NUMB)
C     IF(JJ.GE.99) GO TO 320
C     KK=HJK(NUMB)
C     LL=HJL(NUMB)
C     SINL=SINTOL(NUMB)
C     FOBS=FO(NUMB)
C     SFCALC=SFCA(NUMB)
C     COSAA=COSPHA(NUMB)
C     SINAA=SINPHA(NUMB)
C     DEL=FOBS-SFCALC
      READ(NDATAP,END=320)JJ,KK,LL,SINL,FOBS,SFCALC,COSAA,SINAA,DEL,SSIG
C     IF(EOF(NDATAP).NE.0.0) GO TO 320
      IF(NTAPIN-10)141,6124,141
 141  IF (NTAPSV-100)1004,1004,1005
 1005 FOBS=SQRT(FOBS)
      SFCALC=SQRT(SFCALC)
      DEL=FOBS-SFCALC
      GO TO 1004
 1004 IF (NTAPIN-21)183,184,185
  183 AA=DEL*COSAA
      BB=DEL*SINAA
      GO TO 204
 184  IF (FOBS)183,180,183
  185 IF (NTAPIN-23)186,187,154
  186 AA=FOBS*COSAA
      BB=FOBS*SINAA
      GO TO 204
  187 AA=SFCALC*COSAA
      BB=SFCALC*SINAA
      GO TO 204
 154  IF (NTAPIN-25)155,156,2989
 155  AA=FOBS**2
      BB=0.0
      GO TO 204
 156  AA=SFCALC**2
      BB=0.0
      GO TO 204
 2989 AA=FOBS**2*SFCALC*COSAA
      BB=FOBS**2*SFCALC*SINAA
      GO TO 204
  199 GO TO 201
 6123 GO TO 201
 6124 IF (FOBS)6125,6123,6125
6125  AA=FOBS
      BB=0.0
      GO TO 204
  202 READ (IN,203)JJ,KK,LL,AA,BB
      IF (JJ-999)204,320,320
  204 IF (AA)206,205,206
  205 IF (BB)206,201,206
  206 IF (NREJ)212,210,212
  212 CALL REJECT(JJ,KK,LL,SINL,FOBS,SFCALC,COSAA,SINAA,DEL,NSKIP,AA,BB,
     1 SNLMAX)
      IF (NSKIP)201,210,201
  210 NTOTAL=NTOTAL+1
      IF (JJ)211,207,211
  207 IF (KK)211,208,211
  208 IF (LL)211,209,211
  209 AA=AA/2.
      BB=BB/2.
  211 DO 254 N=1,NEQUIV
  217 J(N)=M1(N,1)*JJ+M2(N,1)*KK+M3(N,1)*LL
      K(N)=M1(N,2)*JJ+M2(N,2)*KK+M3(N,2)*LL
      L(N)=M1(N,3)*JJ+M2(N,3)*KK+M3(N,3)*LL
      NU=J(N)*MTX(N)+K(N)*MTY(N)+L(N)*MTZ(N)
      CALL TRIG(XCOS,XSIN,NU)
      A(N)=AA*XCOS-BB*XSIN
      B(N)=BB*XCOS+AA*XSIN
      IF (ISST(2)-1)224,3012,3012
 3012 FJX=J(N)
      FKX=K(N)
      FLX=L(N)
      WALJ=(U(1,1)*FJX+U(1,2)*FKX+U(1,3)*FLX)
      WALK=(U(2,1)*FJX+U(2,2)*FKX+U(2,3)*FLX)
      WALL=(U(3,1)*FJX+U(3,2)*FKX+U(3,3)*FLX)
      J(N)=WALJ+SIGN(0.5,WALJ)
      K(N)=WALK+SIGN(0.5,WALK)
 3010 L(N)=WALL+SIGN(0.5,WALL)
  224 IF (NORIEN-10)238,225,231
  225 CALL ROTATE(K(N),L(N))
      GO TO 238
  231 CALL ROTATE(J(N),L(N))
  238 IF(J(N)) 239,243,257
  239 J(N)=-J(N)
  240 K(N)=-K(N)
  241 L(N)=-L(N)
      B(N)=-B(N)
      IF(J(N).NE.0) GO TO 257
  243 IF(K(N)) 240,244,257
  244 IF(L(N)) 241,257,257
  257 IF (JMAX-J(N))247,249,249
  247 JMAX=J(N)
  249 IF (KMAX-K(N))250,252,252
  250 KMAX=K(N)
      GO TO 254
  252 IF (KMIN-K(N))254,254,253
  253 KMIN=K(N)
  254 CONTINUE
  260 NEW=NEQUIV-1
      DO 270 NN=1,NEW
      NOW=NN+1
      DO 270 N=NOW,NEQUIV
      IF (A(NN))265,264,265
  264 IF (B(NN))265,270,265
  265 IF (J(NN)-J(N))270,266,270
  266 IF (K(NN)-K(N))270,267,270
  267 IF (L(NN)-L(N))270,268,270
  268 A(N)=0.
      B(N)=0.
  270 CONTINUE
      NTRAN=0
      DO 300 N=1,NEQUIV
      IF (A(N))298,285,298
  285 IF (B(N))298,300,298
  298 AMAX=AMAX+ABS(A(N))
      BMAX=BMAX+ABS(B(N))
      NTRAN=NTRAN+1
      JX(NTRAN)=J(N)
      KX(NTRAN)=K(N)
      LX(NTRAN)=L(N)
      AX(NTRAN)=A(N)
      BX(NTRAN)=B(N)
  300 CONTINUE
      IF((NTRAN+NHOLD).GT.24) GO TO 310
  311 DO 301 N=1,NTRAN
      NHOLDN=NHOLD+N
      JHALF(NHOLDN)=JX(N)
      KHALF(NHOLDN)=KX(N)
      LHALF(NHOLDN)=LX(N)
      AHALF(NHOLDN)=AX(N)
  301 BHALF(NHOLDN)=BX(N)
      NHOLD=NHOLD+NTRAN
      NEXPAN=NEXPAN+NTRAN
      RESID=RESID+FLOAT(NTRAN)*DEL*DEL
      IF(NREF.LE.0) GO TO 201
      NREF=NREF-1
      WRITE(III,305) (JX(N),KX(N),LX(N),AX(N),BX(N),N=1,NTRAN)
      WRITE(III,306)
      GO TO 201
  310 IF(NHOLD.EQ.24) GO TO 397
      NHOLD=NHOLD+1
      DO 339 NN=NHOLD,24
      JHALF(NN)=0
      KHALF(NN)=0
      LHALF(NN)=0
      AHALF(NN)=0.0
  339 BHALF(NN)=0.0
  397 WRITE( IXX) (AHALF(N),BHALF(N),N=1,24),(JHALF(N),KHALF(N),
     Z LHALF(N),N=1,24)
      IF(NHOLD.EQ.0) GO TO 329
      NHOLD=0
      GO TO NXIT,(311,312)
  312 DO 351 N=1,24
  351 JHALF(N)=999
      GO TO 397
33986 WRITE(III,3395)
 3395 FORMAT('0+++WARNING-HIT EOF ON INPUT FILE+++'//)
  320 ASSIGN 312 TO NXIT
      GO TO 310
  329 REWIND IXX
      NCD=(JMAX+1)*(KMAX-KMIN+1)*2
      NEF=(JMAX+1)*NY*2
      IF(NTAPIN.LT.20.OR.NTAPIN.GT.21) GO TO 7009
      RESID=SQRT(2.0*RESID)/VOL
      WRITE(III,7029) RESID
 7029 FORMAT('0SIGMA(RHO) =',F6.3,' E/AEX3 (+2 IF NON CENTRIC)'/)
 7009 IF (NTOTAL)1000,1000,1001
 1000 WRITE (III,2000)
      STOP
 1001 IF (NEXPAN)1002,1002,1003
 1002 WRITE (III,2001)
      STOP
 1003 WRITE (III,352)NTOTAL,NEXPAN
      NEFTBL=NMEMMX-NCD
      WRITE (III,355)NCD
      IF (NEFTBL)4007,4007,3009
 3009 IF (NEFTBL-NEF)360,356,356
  356 WRITE (III,357)NEF
      GO TO 385
  360 NNN=1+(NEF/NEFTBL)
      NY=NEFTBL/(2*(JMAX+1))
      WRITE (III,361)NEFTBL,NNN
 385  RETURN
 4007 WRITE (III,4008)
   10 STOP
      END


      SUBROUTINE ROTATE(A,B)
      Z=A
      A=B
      B=Z
      RETURN
      END


      SUBROUTINE SETU(V,A,B,C,COSA,COSB,COSC)
      DIMENSION DCP(6),U1(3)
      DIMENSION U2(3),U3(3)
      DIMENSION U4(3),U(3,3)
      DIMENSION DM(3,3),RM(3,3)
      DIMENSION V(3,3)
      DCP(1)=A
      DCP(2)=B
      DCP(3)=C
      DCP(4)=COSA
      DCP(5)=COSB
      DCP(6)=COSC
      CALL STOAA(DCP,DM)
      CALL MI3(DM,RM,DETA)
      DO 10 J=1,3
      U1(J)=0.0
 10   U2(J)=0.0
      U1(2)=1.0
      U2(3)=1.0
      CALL NORM(DM,U1)
      CALL NORM(DM,U2)
      CALL CROSS(RM,U1,U2,U3)
      CALL NORM(DM,U3)
      CALL CROSS(RM,U3,U1,U4)
      CALL NORM(DM,U4)
      DO 11 I=1,3
      DO 11 J=1,3
 11   U(I,J)=0.0
      DO 12 J=1,3
      U(1,J)=U3(J)
      U(2,J)=U1(J)
 12   U(3,J)=U4(J)
      CALL MI3(U,V,DETA)
      RETURN
      END


      SUBROUTINE SETUP(VOL)
      DIMENSION TRGTBL(1501),ISST(10)
      DIMENSION WORDS(24)
      COMMON /PEAK/ RMIN,MWG,WGM,WGN,WGO
      COMMON TRGTBL,NORIEN,NXX,NYX,NZ,XMIN,XMAX,YMIN,YMAX,ZMIN
      COMMON ZMAX,DX,DY,DZ,AMAX,BMAX,JMAX,KMAX,KMIN,NMEMMX,NEGINA
      COMMON NCD,NEF,NEFTBL,NX,NY,NCUE,N1,RMAX,NLAST,NFPLT,YFL
      COMMON N11,IN,III,IVXX,IXX,NINPUT,NN3,NDATAP,ISST,A,B,C
      COMMON COSA,COSB,COSC,NSEC,NGRID,NSKIP,CNTRNO,CNTNEG,NEXPAN
      COMMON WORDS
 99   FORMAT (20A4)
 100  FORMAT (1H120A4/)
 2200 FORMAT(10I1,F5.0,E15.6)
 2201 FORMAT('0CONTROL INTEGERS ',10I1,3X'RMIN='F6.0,3X'RMAX='E12.4)
 1000 FORMAT (3F10.4,3F10.6)
 3011 FORMAT('0A=',F8.3,3X'B=',F8.3,3X'C=',F8.3,3X'COSA=',F9.5,3X'COSB='
     A ,F9.5,3X'COSG=',F9.5)
 5081 FORMAT('0UP TO 'I4' POSITIVE PEAKS WILL BE LISTED')
 5012 FORMAT (29H0COMPUTE A FOURIER SUMMATION.)
 5020 FORMAT (19H0PRINT THE RESULTS.)
 5045 FORMAT ('0THERE WILL BE AN ALPHAMERIC FOURIER PLOT')
 5102 FORMAT (34H0TRANSFORM CELL BEFORE RIP AND SIG)
      IVXX=1
      NDATAP=2
      III=61
      IN=60
      IXX=8
      NMEMMX=5900
      REWIND NDATAP
      REWIND IXX
C     CALCULATE COSINE TABLE (TRGTBL) FOR USE IN SIG100
   52 ARG=0
      DARG=.0062831853E0
      DO 15 I=1,1001
      TRGTBL(I)=COS(ARG)
   15 ARG=ARG+DARG
      TRGTBL(1)=1.
      TRGTBL(251)=0.
      TRGTBL(501)=-1.
      TRGTBL(751)=0.
      TRGTBL(1001)=1.
      READ(IN,99,END=24537) (WORDS(I),I=1,20)
C     IF(EOFCKF(IN).EQ.1)STOP
      WRITE (III,100)(WORDS(I),I=1,20)
      READ(IN,2200,END=24537) ISST(1),ISST(2),ISST(3),ISST(6),ISST(5),
     A ISST(4),ISST(7),ISST(8),ISST(9),ISST(10),RMIN,RMAX
C     IF(EOFCKF(IN).EQ.1)STOP
      WRITE(III,2201)         ISST(1),ISST(2),ISST(3),ISST(6),ISST(5),
     A ISST(4),ISST(7),ISST(8),ISST(9),ISST(10),RMIN,RMAX
      READ (IN,1000)A,B,C,COSA,COSB,COSC
      VOL=A*B*C*SQRT(1.0-COSA**2-COSB**2-COSC**2+2.0*COSA*COSB*COSC)
 1613 WRITE (III,3011)A,B,C,COSA,COSB,COSC
 5005 IF (ISST(1))5010,5015,5010
 5010 WRITE (III,5012)
 5015 IF (ISST(3))5019,5053,5019
 5019 WRITE (III,5020)
 5053 IF (ISST(4))5055,5025,5055
 5055 IJA=40*ISST(4)
 5080 WRITE (III,5081)IJA
 5025 IF (ISST(5))5041,1998,5041
 5041 WRITE (III,5045)
 1998 IF (ISST(2))5100,5101,5100
 5100 WRITE (III,5102)
 5101 IF(ISST(7).EQ.0) GO TO 7101
      IJA=4*ISST(7)
      WRITE(III,8101) IJA
 8101 FORMAT('0UP TO 'I4' PEAKS WILL BE PUNCHED')
 7101 IF(ISST(8).EQ.0) GO TO 6200
      IJA=4*ISST(8)
      WRITE(III,9101) IJA
 9101 FORMAT('0UP TO 'I4' NEGATIVE TROUGHS WILL BE LISTED')
 6200 RETURN
24537 STOP
      END


      SUBROUTINE SIG100
Csfy  DIMENSION J(24), K(24), L(24), A(24), B(24), R( 101), TABLE( 5900)
      DIMENSION L(24), A(24), B(24), R( 101), TABLE( 5900)
      DIMENSION WORDS(24)
      INTEGER JH,KH,LH
      DIMENSION JH(24),KH(24),LH(24)
      COMMON TRGTBL(1501), NORIEN,NXMAX,NYMAX,NZ,XMIN
      COMMON  XMAX   , YMIN   , YMAX   , ZMIN   , ZMAX   , DX
      COMMON  DY     , DZ     , AMAX   , BMAX   , JMAX   , KMAX
      COMMON  KMIN   , NMEMMX , NEGINA , NCD    , NEF    , NEFTBL
      COMMON  NX     , NY     , NCUE   , N1     , RMAX   , NLAST
      COMMON  NFPLT  , YFL    , N11    , IN     , III    , IVXX
      COMMON  IXX    , NINPUT , NN3    , NDATAP , ISST(10),AXX
      COMMON  BXX    , C      , COSA   , COSB   , COSC   , NSEC
      COMMON NGRID,NSKIP,CNTRNO,CNTNEG,NEXPAN,WORDS
      COMMON TABLE
  122 FORMAT( 23H0FOURIER (SIG) LEVEL =   I4 )
  441 FORMAT ( 8H0RMAX =   E13.6   ////)
      SEV      = 1000.0
      ZER      =    0.0
      AOM7     =-1000.0
      NUMB=0
      RMAX     =    0.0
      NJ       = JMAX + 1
      NK       = KMAX - KMIN + 1
      NC       = NCD / 2
      Z        = ZMIN
      NZMAX    = ( ( ZMAX - ZMIN ) / DZ + 1.5 )
      NYMAX    = ( ( YMAX - YMIN ) / DY + 1.5 )
      NZ = 1
  100 DO 120 I = 1, NMEMMX
  120 TABLE(I) = 0.0
      WRITE(III,122)NZ
  130 READ(IXX) (A(I),B(I),I=1,24),(JH(I),KH(I),LH(I),I=1,24)
      DO 143 I = 1, 24
      IF(JH(I) .EQ. 999 ) GO TO 175
      LOCC =JH(I) * NK + (KH(I) - KMIN ) + 1
      LOCD = LOCC + NC
      L(I)=LH(I)
      XL = L(I)
      XLZ = XL * Z
      CALL SINCOS ( XLZ, SINLZ, COSLZ )
      TABLE(LOCC) = TABLE(LOCC) + A(I) * COSLZ + B(I) * SINLZ
      TABLE(LOCD) = TABLE(LOCD) - A(I) * SINLZ + B(I) * COSLZ
  143 CONTINUE
      GO TO 130
  175 REWIND IXX
      YI  =  YMIN
      YF  =  YMAX
  204 NY = ( ( ( YF - YI ) / DY ) + 1.5 )
      NEF = ( JMAX + 1 ) * NY * 2
      IF ( (NCD+NEF) .LE. NMEMMX ) GO TO 205
      YF = YF - DY
      GO TO 204
  205 NE = NY * NJ
      XKMIN = KMIN
      LOCEI = NCD + 1
      LOCEF = LOCEI + 2*NE - 1
      DO 206 I = LOCEI, LOCEF
  206 TABLE(I) = 0.0
      LOCC = 1
      LOCD = LOCC + NC
      DO 250 N = 1, NJ
      RKDY = XKMIN * DY
      RKYI = XKMIN * YI
      DO 240 M = 1, NK
      LOCE = LOCEI
      RKY = RKYI
      LOCF = LOCE + NE
      CHY = TABLE(LOCC)
      DHY = TABLE(LOCD)
      DO 230 I = 1, NY
      CALL SINCOS( RKY, SINKY, COSKY )
      TABLE(LOCE) = TABLE(LOCE) + CHY*COSKY + DHY*SINKY
      TABLE(LOCF) = TABLE(LOCF) - CHY*SINKY + DHY*COSKY
      LOCE = LOCE + NJ
      LOCF = LOCE + NE
      RKY = RKY + RKDY
  230 CONTINUE
      RKYI = RKYI + YI
      RKDY = RKDY + DY
      LOCC = LOCC + 1
      LOCD = LOCC + NC
  240 CONTINUE
      LOCEI = LOCEI + 1
  250 CONTINUE
      NX = ( ( XMAX - XMIN ) / DX + 1.5 )
      LOCE = NCD + 1
      LOCF = LOCE + NE
      Y = YI
      WRITE(IVXX) SEV, SEV            ,
     W            NX, NY, XMIN, XMAX, DX, ZMAX, Z, NORIEN, YMIN, YMAX,
     X            DY, AMAX, BMAX, DZ, ZMIN, NZ
      NUMB=NUMB+1
      DO 380 I = 1, NY
      DO 320 N = 1, NX
  320 R(N) = 0.0
      HDX = 0.0
      HXI = 0.0
      DO 360 M = 1, NJ
      E = TABLE( LOCE )
      F = TABLE ( LOCF )
      HX = HXI
      DO 340 N = 1, NX
      CALL SINCOS ( HX, SINHX, COSHX )
      R(N) = R(N) + E * COSHX + F * SINHX
      HX = HX + HDX
  340 CONTINUE
      HDX = HDX + DX
      HXI= HXI+ XMIN
      LOCE = LOCE + 1
      LOCF = LOCE + NE
  360 CONTINUE
      DO 370 N = 1, NX
      IF( R(N) .GT. RMAX ) RMAX = R(N)
  370 CONTINUE
      WRITE ( IVXX ) ZER, ZER       ,
     1             Y, ( R(N), N = 1, NX )
      NUMB=NUMB+1
      Y = Y + DY
  380 CONTINUE
      IF ( YF .GE. YMAX ) GO TO 400
      YI = YF
      YF = YMAX
      GO TO 204
  400 Z = Z + DZ
      NZ = NZ + 1
      IF ( NZ .LE. NZMAX ) GO TO 100
      WRITE ( IVXX ) AOM7, RMAX       ,
     W            NX, NY, XMIN, XMAX, DX, ZMAX, Z, NORIEN, YMIN, YMAX,
     X            DY, AMAX, BMAX, DZ, ZMIN, NZ
      NUMB=NUMB+1
      WRITE(III,45441) NUMB
45441 FORMAT(      '0NO OF BLOCKS USED FOR STORING SUMMATION ANSWERS IS'
     1,I6//)
      WRITE(III, 441) RMAX
      REWIND IVXX
      RETURN
      END


      SUBROUTINE SINCOS ( ARG, SINARG, COSARG )
      COMMON TRGTBL(1501)
    5 IF ( ARG .GE. 0.0 ) GO TO 10
      ARG = ARG + 1.0
      GO TO 5
   10 NARG1 = ( ARG * 1000.0 + 1.5 )
   11 IF ( NARG1 .LE. 1001 ) GO TO 12
      NARG1 = NARG1 - 1000
      GO TO 11
   12 NARG2 = NARG1 + 250
   15 IF ( NARG2 .GT. 1001 ) NARG2 = NARG2 - 1000
   20 SINARG =-TRGTBL ( NARG2 )
      COSARG = TRGTBL ( NARG1 )
      RETURN
      END


      SUBROUTINE STOAA(A,AA)
      DIMENSIONA(6),AA(3,3)
      AA(1,1)=A(1)*A(1)
      AA(2,2)=A(2)*A(2)
      AA(3,3)=A(3)*A(3)
      AA(1,2)=A(1)*A(2)*A(6)
      AA(2,1)=AA(1,2)
      AA(1,3)=A(1)*A(3)*A(5)
      AA(3,1)=AA(1,3)
      AA(2,3)=A(2)*A(3)*A(4)
      AA(3,2)=AA(2,3)
      RETURN
      END


      SUBROUTINE TRANS(U,DETU)
      DIMENSION STUFF(1505)
      DIMENSION ULIMIT(9),EXTRA(19)
      DIMENSION ISST(10),ACELL(6)
      DIMENSION ALSO(562),PLANE(3)
      DIMENSION CELL(3),U(3,3)
      DIMENSION G(3,3),GNEW(3,3)
      DIMENSION UINV(3,3)
      COMMON STUFF,ULIMIT,EXTRA,IN,III,IVXX,IXX,NN5,NN6,NDATAP
      COMMON ISST,ACELL,ALSO,G
   31 FORMAT (9F8.4)
 4    FORMAT ('1TRANSFORMATIONMATRIXREADFROMCARDS....A2=UA1')
    7 FORMAT (7F10.5)
    8 FORMAT (81H1TRANSFORMATION MATRIX TO OBTAIN ORTHOGONAL SYSTEM DEFI
     1NED BY THE FOLLOWING PLANE)
    9 FORMAT (1H010XF10.5,4H X +F10.5,4H Y +F10.5,4H Z =F10.5)
   10 FORMAT (9F8.4)
   11 FORMAT (10H0NEW AXIS ,I2,2X,9HRUNS FROM,F10.5,2X,2HTO,F10.5,2X,
     1 11HIN STEPS OF,F10.5,2X,10HANGSTROMS.)
   17 FORMAT (23H0NEW CELL CONSTANTS ARE6F12.5)
   19 FORMAT (1H0,10X,4HNEWA,I1,1H=,F12.6,5H A1 +,F12.6,5H A2 +,F12.6,
     1 5H A3  )
 1013 FORMAT (35H0THE DETERMINANT OF THE U MATRIX ISF15.6,1H$)
 1001 FORMAT (38H0THE MATRIX INVERSE TO U IS AS FOLLOWS)
 1012 FORMAT (1H0,10X,3F15.6)
      CALL STOAA(ACELL,G)
      IF (ISST(2)-1)1,2,3
    1 RETURN
    2 READ (IN,31)((U(I,J),J=1,3),I=1,3)
C     TRANSFORMATION MATRIX IF IT IS TO BE READ
      WRITE (III,4)
      GO TO 6
    3 READ (IN,7)(PLANE(I),I=1,3),DPLANE,(CELL(J),J=1,3)
C     EQUATION OF PLANE.....SUM(PLANE(I)*X(I))=DPLANE
C     MINIMUM DIMENSIONS OF TRANSFORMED CELL
      WRITE (III,8)
      WRITE (III,9)(PLANE(I),I=1,3),DPLANE
      CALL GP1(PLANE,DPLANE,G,U)
    6 READ (IN,10)(ULIMIT(I),I=1,9)
      IF (ISST(2)-1)1,66,666
  666 DO 12 I=1,3
      WRITE (III,11)I,ULIMIT(2*I-1),ULIMIT(2*I),ULIMIT(I+6)
      ACELL(I)=AMAX1(ABS(ULIMIT(2*I-1))*2.0,ABS(ULIMIT(2*I))*2.0,10.0,
     1CELL(I))
      DO 20 J=1,3
   20 U(I,J)=U(I,J)*ACELL(I)
      ULIMIT(2*I-1)=ULIMIT(2*I-1)/ACELL(I)
      ULIMIT(2*I)=ULIMIT(2*I)/ACELL(I)
      ULIMIT(I+6)=ULIMIT(I+6)/ACELL(I)
   12 ACELL(I+3)=0.0
      GO TO 14
   66 DO 15 I=1,3
      DO 15 J=1,3
      GNEW(I,J)=0.0
      DO 15 K=1,3
      DO 15 L=1,3
   15 GNEW(I,J)=GNEW(I,J)+U(I,K)*U(J,L)*G(K,L)
      DO 16 I=1,3
   16 ACELL(I)=SQRT(GNEW(I,I))
      ACELL(4)=GNEW(1,2)/(ACELL(1)*ACELL(2))
      ACELL(5)=GNEW(1,3)/(ACELL(1)*ACELL(3))
      ACELL(6)=GNEW(2,3)/(ACELL(2)*ACELL(3))
   14 WRITE (III,17)(ACELL(I),I=1,6)
      DO 18 I=1,3
   18 WRITE (III,19)I,(U(I,J),J=1,3)
      CALL MI3(U,UINV,DETU)
      WRITE (III,1013)DETU
      WRITE (III,1001)
      WRITE (III,1012)((UINV(I,J),J=1,3),I=1,3)
      RETURN
      END


      SUBROUTINE TRIG(COS,SIN,NU)
 1170 IF (NU)1173,1201,1171
 1171 NU=NU-12
      GO TO 1170
 1173 NU=NU+12
      IF (NU)1173,1201,1175
 1175 IF (NU-6)1176,1204,1177
 1176 IF (NU-3)1209,1206,1212
 1177 IF (NU-9)1219,1214,1216
 1201 COS=1.
 1202 SIN=0.
      GO TO 1218
 1204 COS=-1.
      GO TO 1202
 1206 SIN=1.
 1207 COS=0.
      GO TO 1218
 1209 COS=.5
 1210 SIN=.866025
      GO TO 1218
 1212 COS=-.5
      GO TO 1210
 1214 SIN=-1.
      GO TO 1207
 1216 COS=.5
 1217 SIN=-.866025
 1218 RETURN
 1219 COS=-.5
      GO TO 1217
      END


      SUBROUTINE VM(A,B,C)
      DIMENSION A(3),C(3)
      DIMENSION B(3,3)
      DO 2 J=1,3
      C(J)=0
      DO 2 I=1,3
    2 C(J)=C(J)+A(I)*B(I,J)
      RETURN
      END


      FUNCTION VMV(A,B,C)
      DIMENSION A(3),B(3,3),C(3)
      VMV=0.0
      DO 1 I=1,3
      DO 1 J=1,3
    1 VMV=VMV+A(I)*B(I,J)*C(J)
      RETURN
      END


      SUBROUTINE REJECT(J,K,L,A,B,C,D,E,F,N,G,H,SNL)
      N=0
      IF(J.EQ.0.AND.K.EQ.0.AND.L.EQ.0) N=1
      IF(A.GT.SNL) N=1
      RETURN
      END

