      SUBROUTINE OPTP (KK,N1,LPIA,LMAXIN,LBORN,U,F,FX)
C *** SEE R H LANDAU, PROGRAM LPOTT, 1981
      IMPLICIT  REAL*8   (A-H,O-Z)
      REAL*8 IMTIJ,MPI,MN,K,KP,KK,IMVABS
	dimension Voptp1(17,17,4,2),Voptp2(17,17,4,2)
  	common /twob/ t2body(17,17,2,14)
      DIMENSION RHO(4,50), BI(100), VCOULL(30), REVABS(30,2), IMVABS(30,
     12)
      DIMENSION NIFTY(20), BNUC(20), BNUCF(20), BN(16), BNF(16), RETIJ
     1(14), IMTIJ(14), KK(50), U(4356,2), F(1089,4), FX(1)
      COMMON /SEC2/ BNUC,BN,BNUCF,BNF,RETIJ,IMTIJ,HBARC,PI,MPI,MN,NZ,NES
     1,NWAVES,NIFTY,NA
      COMMON /SEC4/ ACHP,ACMP,WSP,ACHN,ACMN,WSN
      COMMON /SEC5/ DRHO(4,50),RHO,XKP,XK,IMAX/SEC6/RCOUL,RCUT
	common /nocheck/ nocheck
	data nocheckfirst/1/
      DATA VCOULL,REVABS,IMVABS/150*0.0E+00/
C
C     KK= GRID POINT VECTOR,  N1= SIZE OF MATRIX,KK AND U SIZE SET BY CA
C     CALCULATE PARTIAL WAVE DECOMP OF RHO- THE FORM FACTOR
C     RHO(NFF,I) IS THE I-1 TH PATIAL WAVE AMPLTD OF THE FORM FACTOR
C     FOR A VALUE OF K A ND KP. NFF DETERMINES WHICH F.F. IS USED
C      NFF=1 PROTON MATTER     2 NEUTRON MATTER
C      3 PROTONE SPIN  4 NEUTRON SPIN
C     DO ALL LVALUES ON 1STCALL,THEN STORE ON FILE7, FT IS FOR TEMP
C
      IF (LPIA.NE.1) GO TO 250
C     1ST CALL TO OPTP, CALCULATE ALL U S
      WRITE(6,410) KK(N1)
      K = -200.
      CALL TNCM (K)
      N2 = 2*N1
      NN = NA-NZ
C      MAKE SURE LMAX=1 CASE DOES NOT RUN INTO INDEF ON INTIAL RUNS
      LDMAX = LMAXIN
      IF (LDMAX.EQ.1) LDMAX = 2
      RHL3 = (ACHP/HBARC)**2/2.
      RHL4 = (ACHN/HBARC)**2/2.
C     ADDEDED TERMS FOR C13 SPIN...MAY WANT TO CHANGE SIZE  HERE
      RHL5 = (WSN/HBARC)**2/2.
      AS = 1./6.
      AP = (NZ-2.)/6./NZ*(ACMP/ACHP)**2
      IF (NZ.LE.2) AP = 0.
      AN = (NN-2.)/6./NN*(ACMN/ACHN)**2
      IF (NN.LE.2) AN = 0.
C     SPECIAL C13 , ADD IN P1/2 N AND RESCALE
      IF (NA.EQ.13) AN = (NN-3.)/6./(NN-1.)*(ACMN/ACHN)**2
C     DO LOOP OVER ALL GRID POINTS
C     SIZE OF U MUST BE 4*(NGP+1)**2,2
      REWIND 7
      IF (NIFTY(6).EQ.6) NIFTY(17) = 1
      AOVERA = (NA-1.)/NA
      IF (NIFTY(17).EQ.1) AOVERA = 1.
      IF (NA.EQ.2) AOVERA = 1.
      IF ((NIFTY(6).EQ.4).OR.(NIFTY(6).EQ.6)) AOVERA = 0.
      XHL = AOVERA*2.*PI*PI
      APARAM = 1.
      I = NIFTY(1)+1
      GO TO (30,10,20,20,20,30,20,30,20,20), I
 10   APARAM = -1.
      GO TO 30
 20   APARAM = 0.
 30   APARAM = (1.5*NZ*APARAM/RCOUL)*HBARC/137.036
      XHL1 = RCOUL/HBARC
      XHL2 = RCUT/HBARC
      NIF = NIFTY(10)
      DO 240 I1=1,N1
         DO 240 I2=1,I1
            K = KK(I2)
            KP = KK(I1)
	 	jj1 = I1
		jj2 = I2
            LMAX = LDMAX+5
            IF (NIFTY(16).NE.0) GO TO 90
C      FO SHELL MODEL NUCLESU(C,0) USE ANALYTIC EXPRESSION WITH IL*S
C     DIFFERENT NEUTRON AND PROTONE DENSITEIES POSSIBLE
            ZP = K*KP*RHL3
            XP = (K*K+KP*KP)*RHL3/2.
            CALL BESLIA (ZP,BI)
            B = 0.
            IF ((ZP-XP).GT.-150.) B = EXP(ZP-XP)
            DO 40 I=1,LMAX
               RHO(3,I) = 0.
               RHO(1,I) = B*(2.*I-1)*((1.-4.*AP*(XP-I+1))*BI(I)+4.*ZP*AP
     1         *BI(I+1))
 40         CONTINUE
            IF (NA.EQ.13) GO TO 60
            IF ((ACHP.NE.ACHN).OR.(ACMP.NE.ACMN)) GO TO 60
C     EQUAL N AND P DENSITIES
            DO 50 I=1,LMAX
               RHO(4,I) = 0.
               RHO(2,I) = RHO(1,I)
 50         CONTINUE
            NFF = 1
            GO TO 110
C     UNEQUAL N AND P DENSITES IN H.O. MODEL
 60         NFF = 2
            ZN = K*KP*RHL4
            XN = (K*K+KP*KP)*RHL4/2.
            CALL BESLIA (ZN,BI)
            B = 0.
            IF ((ZN-XN).GT.-150.) B = EXP(ZN-XN)
            DO 70 I=1,LMAX
               RHO(2,I) = B*(2.*I-1)*((1.-4.*AN*(XN-I+1))*BI(I)+4.*ZN*AN
     1         *BI(I+1))
               RHO(4,I) = 0.
 70         CONTINUE
            IF (NA.NE.13) GO TO 110
C     SPIN FOR C13
            NFF = 4
            ZS = K*KP*RHL5
            XS = (K*K+KP*KP)*RHL5/2.
            CALL BESLIA (ZS,BI)
            B = 0.
            IF ((ZS-XS).GT.-150.) B = EXP(ZS-XS)
C     NORMALIZE TO 1/NN SO CAN MULT BY NN AND HAVE 1 NEUTRON
            B = B/NN
            DO 80 I=1,LMAX
               RHO(4,I) = B*(2.*I-1)*((1.-4.*AS*(XS-I+1))*BI(I)+4.*ZS*AS
     1         *BI(I+1))
               RHO(2,I) = ((NN-1.)/NN)*RHO(2,I)+RHO(4,I)
 80         CONTINUE
            IF ((I1.EQ.I2).AND.(I1.EQ.N1)) WRITE (6,350) K,(RHO(4,I),I=1
     1      ,LMAX)
            GO TO 110
C      FOR WOOD SAXON TYPR(OR ARBITRARY) NUCLERA N DENSITY,CALC RHO NUME
 90         XKP = KP
            XK = K
            IMAX = LDMAX+5
            LMAX = IMAX
C     DETEERMINE NUMBER OF F.F.*S NEEDED,IF NIFTY(6)=3, NEED SPIN IN
            NFF = 2
            IF ((ACHP.EQ.ACHN).AND.(ACMP.EQ.ACMN).AND.(WSP.EQ.WSN)) NFF
     1      = 1
            IF (NIFTY(6).EQ.3) NFF = 4
C      FOR RHO**2 IN CA
            IF (NIFTY(16).EQ.2) NFF = 4
            CALL RHOKFF (NFF)
            IF (NFF.NE.1) GO TO 110
C     NFF=1, SEET N AND P MATTER DENSITES=, ALL SPIN DENS=0
            DO 100 I=1,LMAX
               RHO(2,I) = RHO(1,I)
               RHO(3,I) = 0.
               RHO(4,I) = 0.
 100        CONTINUE
 110        CONTINUE
            IF ((NIFTY(10).NE.3).AND.(NIFTY(10).NE.4)) GO TO 120
C     INCLUDE CUTOFF COULOMB POTENTIAL IN MOMENTUM SPACE
            AK = K
            AKP = KP
            CALL VCOUL (AK,AKP,APARAM,XHL1,XHL2,VCOULL,LDMAX,NIF)
 120  IF (NIFTY(7) .LT. 3 .AND. NIFTY(7).NE.5) GO TO 130
C      INCLUDE ABSORPTION ON 2 NUCLEONS
            CALL VABS (K,KP,KK(N1),REVABS,IMVABS,LDMAX)
 130        CONTINUE
C     NEW VERSION....OUTTER DO LOOPN OVER KP,K,,INNEROVER ALL L S
            XMN = MN
	    CALL TNUC (jj1,jj2,KP,K,KK(N1),XMN)
c----------------------------------------------------------------------
C     DO LOOP OVER L                                                 OPT
C
            DO 230 LDUM=1,LDMAX
               IF ((LDUM.GT.LBORN).AND.((I1+I2).NE.(2*N1))) GO TO 230
               J1 = LDUM
               URC = 0.
               UIC = 0.
               URS = 0.
               UIS = 0.
               IF (NIFTY(6).EQ.3) RHL1 = J1*(J1-1.)
C     DO OVER PI-N ANGULAR MOMENTUM;S,P,D,F,G,
               DO 200 J2=1,5
                  J4 = IABS(J1-J2)+1
                  J5 = J1+J2-1
                  IF (NIFTY(6).EQ.3) RHL2 = J2*(J2-1.)
                  DO 190 J3=J4,J5
                     RHL = CGC2(J3,J2,J1)
C     SPIN INDEPENDT POTENTIAL
                     IF (NIFTY(1).EQ.1.OR.NIFTY(1).EQ.8) GO TO 140
C     PI PLUS SCATTERINF
                     JPR = 4*J2-3
                     JPI = JPR+1
                     JNR = JPR+2
                     JNI = JPR+3
                     GO TO 150
 140                 CONTINUE
C     PI MINUS SCATTERIN, REARRANGE CODE
                     JNR = 4*J2-3
                     JNI = JNR+1
                     JPR = JNR+2
                     JPI = JNR+3
 150                 CONTINUE
                     IF (NIFTY(1).NE.4) GO TO 160
C     PI ZERO SCATTERING,  TAKE AVERAGE OF P AND N AMPLITUDES
                     BNUC(JPR) = (BNUC(JPR)+BNUC(JNR))/2.
                     BNUC(JNR) = BNUC(JPR)
                     BNUC(JPI) = (BNUC(JPI)+BNUC(JNI))/2.
                     BNUC(JNI) = BNUC(JPI)
                     IF (NIFTY(6).NE.3) GO TO 160
                     BNUCF(JPR) = (BNUCF(JPR)+BNUCF(JNR))/2.
                     BNUCF(JNR) = BNUCF(JPR)
                     BNUCF(JPI) = (BNUCF(JPI)+BNUCF(JNI))/2.
                     BNUCF(JNI) = BNUCF(JPI)
 160                 CONTINUE
                     URC = URC+RHL*(NZ*RHO(1,J3)*BNUC(JPR)+NN*RHO(2,J3)*
     1               BNUC(JNR))
                     UIC = UIC+RHL*(NZ*RHO(1,J3)*BNUC(JPI)+NN*RHO(2,J3)*
     1               BNUC(JNI))
                     IF ((NIFTY(3).NE.2).AND.(NIFTY(3).NE.3)) GO TO 180
C      PRINT OUT OFF-SHELL TMATRIX ONTO TAPE 3
                     IF (LDUM.NE.1) GO TO 180
                     IF ((I1.NE.I2).OR.(I1.NE.1)) GO TO 170
C     FIRST TIME THRU, PRINT OUT IDENTIFIERS
                     REWIND 3
                     WRITE (3,360) NZ,NA,(NIFTY(N),N=1,20)
                     WRITE (3,370) ACHP,ACMP,WSP,ACHN,ACMN,WSN,RCOUL,
     1               RCUT
                     NGP = N1-1
                     WRITE (3,360) NGP
                     WRITE (3,380) (KK(N),N=1,N1)
C     PRINT OUT PI N T MATRIX IN PI-A COM  FOR K' GTR K
 170                 IF (I1.LT.I2) GO TO 180
                     WRITE (3,390) I1,I2,BNUC(JPR),BNUC(JPI),BNUC(JNR),
     1               BNUC(JNI)
                     IF (NIFTY(6).EQ.3) WRITE (3,380) BNUCF(JPR),BNUCF
     1               (JPI),BNUCF(JNR),BNUCF(JNI)
 180                 IF (NIFTY(6).NE.3) GO TO 190
C     ---- ---SPIN DEPENDNET POTENTIAL
                     IF (J1.EQ.1) GO TO 190
                     RHL = RHL*(RHL2+RHL1-J3*(J3-1.))/(2.*RHL1)
                     URS = URS+RHL*(NZ*RHO(3,J3)*BNUCF(JPR)+NN*RHO(4,J3)
     1               *BNUCF(JNR))
                     UIS = UIS+RHL*(NZ*RHO(3,J3)*BNUCF(JPI)+NN*RHO(4,J3)
     1               *BNUCF(JNI))
 190              CONTINUE
 200           CONTINUE
C     FORM POTENTIALS FOR STATE OF DEFINITE J
c********* Mr Lu stores optp here and prints it    later *********
c pure unclear part
        if (LDUM .le. 4) then
                Voptp1(I1,I2,LDUM,1) = URC
                Voptp1(I1,I2,LDUM,2) = UIC
        endif

               RHL = XHL/(2.*J1-1.)
               VC = 0.
               IF ((NIFTY(10).NE.3).AND.(NIFTY(10).NE.4)) GO TO 210
               VC = VCOULL(J1)*AOVERA
               IF (NIFTY(6).EQ.6) VC = VCOULL(J1)
 210           IF (NIFTY(7).NE.3.AND.NIFTY(7).NE.4) GO TO 220
C     ADD IN VABS TO CENTRAL POTENTIAL
C     AOVERA SCALING ALREADY DONE FOR VABS IN SUB VABS
               URC = URC+REVABS(J1,1)/RHL
               UIC = UIC+IMVABS(J1,1)/RHL
	if (LDUM .le. 4) then
                Voptp2(I1,I2,LDUM,1) = REVABS(J1,1)/RHL
                Voptp2(I1,I2,LDUM,2) = IMVABS(J1,1)/RHL
        endif

               IF (NIFTY(7).NE.4.OR.NIFTY(6).NE.3) GO TO 220
C     INCLUDE SPIN-DEPENDENT ABSORPTION
               URS = URS+REVABS(J1,2)/RHL
               UIS = UIS+IMVABS(J1,2)/RHL
 220           CONTINUE
C     SET UP POTENTIAL NEEDED FOR THIS L IN INTEGRAL EQTN
               NPOT1 = 2*LDUM
               F(NPOT1-1,1) = RHL*(URC+(J1-1)*URS)+VC
               F(NPOT1,1) = -RHL*(UIC+(J1-1)*UIS)
               NNN = NPOT1-1
                 IF(KP.EQ.KK(N1).AND.KP.EQ.K.AND.LDUM.EQ.1)
     1           WRITE(6,410)KK(N1)
               IF ((KP.EQ.KK(N1)).AND.(LDUM.EQ.2)) WRITE (6,400) K,F(NNN
     1         ,1),VC,REVABS(J1,1),REVABS(J1,2),F(NPOT1,1),IMVABS(J1,1),
     2         IMVABS(J1,2)
               IF (NIFTY(6).NE.3) GO TO 230
               F(NPOT1-1,2) = RHL*(URC-J1*URS)+VC
               F(NPOT1,2) = -RHL*(UIC-J1*UIS)
C     DO LOOP OVER LDUM ENDS
 230        CONTINUE
            NPOT1M = 2*LBORN
            IF (LDMAX.LT.LBORN) NPOT1M = 2*LDMAX
            IF ((I1.EQ.I2).AND.(I1.EQ.N1)) NPOT1M = 2*LDMAX
C      NEW UNFORMATTED READ FOR TIME SAVING
            IF (NIFTY(6).NE.3) WRITE (7) (F(NPOT1,1),NPOT1=1,NPOT1M)
            IF (NIFTY(6).EQ.3) WRITE (7) ((F(NPOT2-1,I),F(NPOT2,I),I=1,2
     1      ),NPOT2=2,NPOT1M,2)
C     DO LOOP OVER KP AND K ENDS
 240  CONTINUE
c ******** Mr Lu prints out optp here ************
	do ll = 1, 4
		write(6,500) ll, (kk(jj),jj=1,16) 
		do ii = 1, 16
		  write(6,501) (Voptp1(ii,jj,ll,1),jj=1,16)
		enddo
		write(6,502)
		do ii = 1, 16
                  write(6,501) (Voptp1(ii,jj,ll,2),jj=1,16)
                enddo
	enddo
        do ll = 1, 4
                write(6,503) ll, (kk(jj),jj=1,16)
                do ii = 1, 16
                  write(6,501) (Voptp2(ii,jj,ll,1),jj=1,16)
                enddo
                write(6,502) 
                do ii = 1, 16
                  write(6,501) (Voptp2(ii,jj,ll,2),jj=1,16)
                enddo
        enddo

500	format(' **** Helleo Mr Lu **** Angular moment = ', i5,
     &	5x, 'and grid points of momentum:'/16e8.2)   
503     format(' **** Absp, Mr Lu **** Angular moment = ', i5,
     &  5x, 'and grid points of momentum:'/16e8.2)
501	format(1x,16e8.2)
502	format(' Imaginary part: ') 
c---------- print off-shell two body potential -----------
	open(12,file='t2b_pot',status='unknown')
	write(12,505)
505	format(/'Mr Lu print off-shell two body potential:')
	write(12,506) (kk(n), n=1,16)
	do ncom = 1,2
	do j = 1,6
	  do jj1 = 1,16
		write(12,506) (t2body(jj1,jj2,ncom,j), jj2 = 1,16)
	  enddo
	enddo
	enddo
506	format(16e15.8)
 250  REWIND 7
C     POTS ALL STORED,READ IN THE VALUES FOR CUURRENT LPIA
C     NPTS =? OF I1,I2 COMBOS(FOLLOWED BY U S FOR ALL L S)
      NPTS = (N1*(1+N1))/2
      IF (LPIA.NE.1) GO TO 260
C     FIRST CALL FOR READ, SETUP PARAMS
C     N.B....THE VALUES OF VL(I,J) FOR ALL L BUT FIXED
C     I,J ARE ON ONE RECORD VIA UNFORMATTED WRITE
C      NEED SKIP TO NEW RECORD TO GET NEW IJ,
      NSPIN = 1
      NPOT = 2
      IF (NIFTY(6).NE.3) GO TO 260
C     SPIN-DEPENDENT CASE, NOW 4 V S ADJACENT FOR EACH L VS 2 IF NF
      NSPIN = 2
      NPOT = 2*NPOT
 260  LDUM = 2*(LPIA-1)*NSPIN
      IF (LPIA.GT.LBORN) GO TO 290
C     DO LOOP OVER NPTS (I.E. RECORDS)
      DO 280 NPT=1,NPTS
         IF (LPIA.NE.1) GO TO 270
         READ (7,end=280) (F(NPT,I),I=1,NPOT)
         GO TO 280
C        SPACE TO L VALUE NEEDED
 270     READ (7,end=280)(RHL,I=1,LDUM),(F(NPT,I),I=1,NPOT)
 280  CONTINUE
      GO TO 310
C     LPIA .GT. LBORN, READ ONLY LAST ON-SHELL POINT
 290  NSKIP = NPTS-1
      DO 300 I=1,NSKIP
         READ (7)
 300  CONTINUE
      READ (7,end=310) (RHL,I=1,LDUM),(F(NPTS,I),I=1,NPOT)
C     SET UP POTENTIAL MATRIX
 310  ISTART = 1
      IF (LPIA.GT.LBORN) ISTART = N1
      DO 340 I1=ISTART,N1
         DO 340 I2=ISTART,I1
            IFLAG = 0
            NPT = ((I1-1)*I1)/2+I2
            NPOT1 = (I2-1)*N2+I1
 320        NPOT2 = NPOT1+N1*(N2+1)
            NPOT3 = NPOT2-N1
            NPOT4 = NPOT1+N1
            U(NPOT1,1) = F(NPT,1)
            U(NPOT2,1) = U(NPOT1,1)
            U(NPOT3,1) = F(NPT,2)
            U(NPOT4,1) = -U(NPOT3,1)
C     SPIN-DEPENDNT POTENTIALS
            IF (NIFTY(6).NE.3) GO TO 330
            U(NPOT1,2) = F(NPT,3)
            U(NPOT2,2) = U(NPOT1,2)
            U(NPOT3,2) = F(NPT,4)
            U(NPOT4,2) = -U(NPOT3,2)
C     SET UP POTNTIAL WITH U(KP,K)=U(K,KP)..SYMMETRIC U
 330        IF ((I1.EQ.I2).OR.(IFLAG.EQ.1)) GO TO 340
            NPOT1 = (I1-1)*N2+I2
            IFLAG = 1
            GO TO 320
 340  CONTINUE
      RETURN
C
 350  FORMAT (30H C13 NEUTRON SPIN FF FOR K=KP=,E15.5/10(10X,6E15.6/))
 360  FORMAT (2I5,2X,I1,I2,8I1,10I3,F10.3,E13.6)
 370  FORMAT (8F10.4)
 380  FORMAT (5E16.7)
 390  FORMAT (2I5/5E16.7)
 400  FORMAT (3X,3E19.4,2E11.3,E19.4,2E11.3)
 410  FORMAT(28H0OPT POTEN VALUES FOR KP=K0=,F10.3,7HAND L=1/20X,1HK,
     1 13X,3HREV,16X,2HVC,15H RABS(N0 FLIP) ,11H RABS(FLIP),3HIMV,
     217X,27H IMVABS(NOFLIP),IMABS(FLIP)/)
      END
