C                              PROGRAM EXT_J
C**************************************************************************
C PROGRAMS WHICH EXTRACT SCALAR COUPLINGS FROM MEASUREMENT OF SEPARATION
C OF EXTREMA IN DISPERSIVE AND ABSORPTIVE PLOTS OF ROWS THROUGH CROSS PEAKS 
C BY ITERATION METHOD.
C WRITTEN BY YANGMEE KIM ( OCT. 25 1988)
C REFERENCE: J. M. R 84, 9-13 (1989)        
C**************************************************************************
C   
      CHARACTER*30 OUTPUTFILE,INPUTFILE,TITLE
      REAL*8 X(999),Y(999),Z,COUJ(1000)
&           ,WO(1000),DATAJ(1000),DATAW(1000),SN(1000),SW
      INTEGER NR(1000)
      COUPLE(Z,AA,BB,CC)=Z*Z*Z+AA*Z*Z+BB*Z+CC
C
      write(*,*) '****************************************************'
      write(*,*) '* THIS PROGRAM CALCULATES COUPLING CONSTANT FROM   *' 
      write(*,*) '* ABSORPTIVE AND DISPERSIVE PEAK-TO-PEAK SEPARATION*'
      write(*,*) '****************************************************' 
C
C INPUT DATA
C
      write(*,"('Enter the name of the input file:',$)")
      read (*,2) INPUTFILE
 2    FORMAT(A30)
      OPEN (UNIT = 4, FILE = INPUTFILE)
      READ(4,11) TITLE
 11   FORMAT(A20)
            NN=0
 50         NN=NN+1
            READ(4,15,ERR=60)NR(NN),X(NN),Y(NN),SN(NN)
            GO TO 50
 60   CONTINUE
      NN=NN-1
 15   FORMAT(I10,3F10.3)
C******************************************************************
C INPUT FILE
C     NR IS RESIDUE NUMBER
C     X IS ABSORPTIVE PEAK TO PEAK SEPARATION (IN DATA POINT IN 4K)
C     Y IS DISPERSIVE PEAK TO PEAK SEPARATION (IN DATA POINT IN 4K)
C     SN IS SIGNAL TO NOISE RATIO.
C     SW IS THE SWEEP WIDTH OF YOUR SPECTRA.
C     !!!!ZERO FILL YOUR DATA TO 4K(4096)!!!!
C******************************************************************
C
C LET'S CALCULATE COUPLING CONSTANT.
      write(*,"('What is the sweep width of your data?',$)")
      read(*,40) SW
      write(*,*) 'How do you want to calculate J?'
      write(*,*) '1.  Iteration method'
      write(*,*) '2.  Analytical solution'
      write(*,"('-->',$)")
      read(*,40) JI
40    format(F8.1)
      IF(JI.EQ.2) GO TO 1001
      DO 561, I=1,NN
         DATAJ(I)=X(I)*SW/4096.0
         DATAW(I)=Y(I)*SW/4096.0
         Y2=DATAW(I)*DATAW(I)
         X2=DATAJ(I)*DATAJ(I)
         X4=X2*X2
         X6=X2*X4
         Y4=Y2*Y2
         Y6=Y2*Y4
         Y8=Y6*Y2
         A1=-Y2
         B1=-9./4.*X4+1.5*X2*Y2+0.75*Y4
         C1=81./64.*X6-9./16.*X4*Y2-21./32.*X2*Y4-1./16.*Y6
&           +Y8/64./X2
         XX1=4.0
         XX2=400.0
         YY1=COUPLE(XX1,A1,B1,C1)
         YY2=COUPLE(XX2,A1,B1,C1)
 52      IF (YY1.EQ.YY2) GO TO 561
         XR=(XX1*YY2-XX2*YY1)/(YY2-YY1)
         YR=COUPLE(XR,A1,B1,C1)
         DX=ABS(1-XX2/XR)
         IF(DX.LT.1.0E-07) GO TO 54
         XX1=XX2
         YY1=YY2
         XX2=XR
         YY2=YR 
         GO TO 52
  54     IF(XR.LT.0.0)XR=0.0
         COUJ(I)=SQRT(XR)
         OMEGA2=9./4.*X4-0.25*Y4-1.5*XR*X2+0.5*XR*Y2-XR*XR
         OMEGA1=XR-1.5*X2-0.5*Y2
         OMEGA=OMEGA2/OMEGA1
         IF (OMEGA.LT.0.0) GO TO 561
         WO(I)=SQRT(OMEGA)
  561  continue
       go to 1110
 1001  DO 1050, I=1,NN
         DATAJ(I)=X(I)*SW/4096.0
         DATAW(I)=Y(I)*SW/4096.0
         Y2=DATAW(I)*DATAW(I)
         X2=DATAJ(I)*DATAJ(I)
         X4=X2*X2
         X6=X2*X4
         Y4=Y2*Y2
         Y6=Y2*Y4
         Y8=Y6*Y2
         Y10=Y8*Y2
         Y12=Y10*Y2
         Y14=Y12*Y2
         Y16=Y14*Y2
         X8=X6*X2
         X10=X8*X2
         X12=X10*X2   
         A1=-Y2
         B1=-9./4.*X4+1.5*X2*Y2+0.75*Y4
         C1=81./64.*X6-9./16.*X4*Y2-21./32.*X2*Y4-1./16.*Y6
&           +Y8/64./X2
         P=B1-A1*A1/3.
         Q=C1-A1*B1/3.+2.*A1*A1*A1/27.
         SQ=Q*Q/4.+P*P*P/27.
         SQ2=SQRT(SQ)
         AAA=-Q/2.+SQ2
         BBB=-Q/2.-SQ2
         AQ=AAA**(0.33333333)
         BQ=-((-BBB)**(0.33333333))
         ABQ=AQ+BQ
         COUT=ABQ-A1/3.
         COUJ(I)=SQRT(COUT)
         XR=COUT
         OMEGA2=9./4.*X4-0.25*Y4-1.5*XR*X2+0.5*XR*Y2-XR*XR
         OMEGA1=XR-1.5*X2-0.5*Y2
         OMEGA=OMEGA2/OMEGA1
         IF(OMEGA.LT.0.0) GO TO 1050
         WO(I)=SQRT(OMEGA)
 1050    continue 
C     COUJ IS COUPLING CONSTANT IN HZ.
C     WO IS LINE WIDTH IN HZ.
C
C OUTPUT RESULTS
c
 1110 write(*,"(' Enter the name of the the outputfile:',$)") 
      read(*,2) OUTPUTFILE  
      OPEN (UNIT = 3, FILE = OUTPUTFILE)  
      write(3,140)
 140  format('**************************************************',
&              '*****************************')
      WRITE (3,141)
 141  FORMAT(/,'               J COUPLING CONSTANT AND LINE WIDTH',/)
      write(3,143)
 143  format('**************************************************',
&              '*****************************')
      write(3,145)
 145  format(/,'    RES   ',
&              '     ABSORP  ' ,
&              ' DISPER **',' AB IN HZ ','DS IN HZ ',
&              '    J   ','    w     ',' S/N ratio'
&              //,'************************************************',
&               '*******************************',/) 
       DO 302, I = 1,NN
      WRITE(3,144)NR(I),X(I),Y(I),DATAJ(I),DATAW(I)
&                 ,COUJ(I),WO(I),SN(I)
 302  continue 
 144  FORMAT(I7,3X,2f10.1,2F10.2,3F9.1/)
      CLOSE (UNIT = 4)        
      CLOSE (UNIT = 3)
 100  STOP
      END