cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Program
c Microscopic Global Optical Potential for heavy-ion reaction
c (mgop.f)
c
c  2011/03/02   version 1.00
c  2011/05/17   version 1.01
c  2011/08/18   version 1.02
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit real*8(a-h,o-z)
      character*40 input,output
      dimension amp(20),rgauss(10)
      parameter(mp=100)
      dimension X(0:MP),X2(0:MP),yamp(0:mp,8,20)
c input
      open(11,file='mgop.in')
      read(11,'(a40)')input
      read(11,'(a40)')output
      read(11,'(i6,f10.3)')ie,ee
      read(11,'(2f6.2)')amtag,ztag
      read(11,'(2f6.2)')rmax,dr
      close(11)
c read input file
      open(11,file=input)
      read(11,*)int,num,ampro,zpro
      do 3001 m=1,8
      read(11,*)tmass
      if(m.eq.1) tagmin=tmass
      if(m.eq.8) tagmax=tmass
      x2(m-1)=tmass**(1./3.)
      do 3002 i=0,num-1
      read(11,*)x(i),
     & yamp(i,m, 1),yamp(i,m, 2),yamp(i,m, 3),yamp(i,m, 4),yamp(i,m, 5),
     & yamp(i,m, 6),yamp(i,m, 7),yamp(i,m, 8),yamp(i,m, 9),yamp(i,m,10),
     & yamp(i,m,11),yamp(i,m,12),yamp(i,m,13),yamp(i,m,14),yamp(i,m,15),
     & yamp(i,m,16),yamp(i,m,17),yamp(i,m,18),yamp(i,m,19),yamp(i,m,20)
 3002 continue
 3001 continue
      close(11)
c check input parameter
      if(ie.ne.1.and.ie.ne.2) then
      write(*,*)'Error! ie=',ie
      stop
      end if
      if(ie.eq.1)ene=ee
      if(ie.eq.2)ene=ee/ampro
      if(int.eq.1.and.(ene.lt.50.d0.or.ene.gt.400.d0)) then
      write(*,*)'Error! int=',int,'ee=',ene,'MeV/u'
      stop
      end if
      if(int.eq.2.and.(ene.lt.200.d0.or.ene.gt.400.d0)) then
      write(*,*)'Error! int=',int,'ee=',ene,'MeV/u'
      stop
      end if
      if(idnint(amtag).lt.idnint(tagmin).or.
     &   idnint(amtag).gt.idnint(tagmax)) then
      write(*,*)'Error! amtag=',amtag,'(tagmax,tagmin)=',tagmax,tagmin
      stop
      end if
c search parameter
      call search(num,ene,amtag,ztag,amp,mp,x,x2,yamp)
      rgauss(10)=0.45d0*(ampro**(1./3.)+amtag**(1./3.)+1.d0)
      rgauss( 9)=rgauss(10)*17./18.
      rgauss( 8)=rgauss(10)*16./18.
      rgauss( 7)=rgauss(10)*15./18.
      rgauss( 6)=rgauss(10)*14./18.
      rgauss( 5)=rgauss(10)*13./18.
      rgauss( 4)=rgauss(10)*12./18.
      rgauss( 3)=rgauss(10)*11./18.
      rgauss( 2)=rgauss(10)*10./18.
      rgauss( 1)=rgauss(10)* 9./18.
c make potential & write
 1000 format('# Microscopic global optical potential with CEG07b')
 1001 format('# Microscopic global optical potential with CEG07a')
 1002 format('# Incident energy (E/A) = ',f12.5,'(MeV/u)')
 1003 format('# projectile    (A,Z,N) = (',f5.0,',',f5.0,',',f5.0,')')
 1004 format('# target        (A,Z,N) = (',f5.0,',',f5.0,',',f5.0,')')
 1005 format('# ')
 1006 format('# R(fm)    V_real(MeV)     W_imag(MeV)')
 1007 format(f8.4,1p2d16.6)
      open(11,file=output)
      if(int.eq.1) write(1,1000)
      if(int.eq.2) write(1,1001)
      write(11,1002)ene
      write(11,1003)ampro,zpro,ampro-zpro
      write(11,1004)amtag,ztag,amtag-ztag
      write(11,1005)
      write(11,1006)
      nrrmax=rmax/dr+1.0001
      do 100 i=1,nrrmax
      r=(i-1)*dr
      v=0.d0
      w=0.d0
      do 200 n=1,10
      v=v+amp(n   )*exp(-r**2/rgauss(n)**2)
      w=w+amp(n+10)*exp(-r**2/rgauss(n)**2)
 200  continue
      write(11,1007)r,v,w
 100  continue
      close(11)
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine search(num,ene,amtag,ztag,amp,mp,x,x2,yamp)
      implicit real*8(a-h,o-z)
      dimension amp(20)
      dimension X(0:MP),Y(0:MP),X2(0:MP),Y2(0:MP)
      dimension yamp(0:mp,8,20),y2amp(0:mp,20)
      a=amtag**(1./3.)
c energy
      do 1001 m=1,8
      do 1003 n=1,20
      do 1004 i=0,num-1
      y(i)=yamp(i,m,n)
 1004 continue
      call spline(num,ene,yy,x,y,mp)
      y2amp(m-1,n)=yy
 1003 continue
 1001 continue
c target
      do 2001 n=1,20
      do 2002 m=1,8
      y2(m-1)=y2amp(m-1,n)
 2002 continue
      call spline(8,a,yy,x2,y2,mp)
      amp(n)=yy
 2001 continue
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine spline(N2,xx,yy,x,y,mp)
      implicit real*8(a-h,o-z)
c
c     N2:number of points(input)
c     xx:your wish point (input)
c     yy:your wish point (output)
c     x :x component(input)
c     y :y component(input)
c     mp:dimension of x and y components(input)
C
C     SPLINE METHOD
C
      DIMENSION X(0:MP),Y(0:MP),G(0:MP),H(0:MP)
C
      N=N2-2
      CALL SPL21(X,Y,G,H,N,MP)
C     DO 40 I=0,N+1
C  40 WRITE(6,*) X(I),Y(I),G(I),H(I)
      YY=SPL22(G,H,X,Y,XX,N,MP)
      return
      END
C
C     CALCULATION OF F(X) WITH CUBIC SPLINE METHOD
C
      FUNCTION SPL22(G,H,X,Y,XX,N,MP)
      implicit real*8(a-h,o-z)
      DIMENSION X(0:MP),Y(0:MP),G(0:MP),H(0:MP)
      DO 10 I=1,N+1
      IF(X(I-1).LE.XX.AND.XX.LE.X(I)) GO TO 20
   10 CONTINUE
   20 XI=XX-X(I-1)
      XL=X(I)-X(I-1)
      SPL22=Y(I-1)+XI*(G(I-1)
     1     +0.5*XI*(H(I-1)+XI*(H(I)-H(I-1))/(3.0*XL)))
      RETURN
      END
C
C     CALCULATION OF THE COEFFICIENTS
C
      SUBROUTINE SPL21(X,Y,G,H,N,MP)
      implicit real*8(a-h,o-z)
      PARAMETER (NP=100)
      DIMENSION A(NP,NP+1),X(0:MP),Y(0:MP),G(0:MP),H(0:MP)
      N1=N+1
      DO 10 I=1,N
      DO 10 J=1,N
   10 A(J,I)=0.0
      DO 11 I=1,N
      XL=X(I)-X(I-1)
      XR=X(I+1)-X(I)
   11 A(I,N1)=6.0*(Y(I-1)/XL-(XL+XR)/XL/XR*Y(I)+Y(I+1)/XR)
      A(1,1)=2.0*(X(2)-X(0))
      A(1,2)=X(2)-X(1)
      A(N,N)=2.0*(X(N1)-X(N-1))
      A(N,N-1)=X(N)-X(N-1)
      DO 12 I=2,N-1
      A(I,I-1)=X(I)-X(I-1)
      A(I,I)=2.0*(X(I+1)-X(I-1))
   12 A(I,I+1)=X(I+1)-X(I)
      CALL GASJO(A,DET,N,NP)
      DO 13 I=1,N
   13 H(I)=A(I,N1)
      H(0)=0.0
      H(N1)=0.0
      XL=X(1)-X(0)
      G(0)=(Y(1)-Y(0))/XL-XL*H(1)/6.0
      DO 14 I=1,N1
   14 G(I)=G(I-1)+0.5*(H(I)+H(I-1))*(X(I)-X(I-1))
      RETURN
      END
C
C     GAUSS-JORDAN METHOD
C
      SUBROUTINE GASJO(A,DET,N,NP)
      implicit real*8(a-h,o-z)
C
      PARAMETER (MP=100)
      DIMENSION A(NP,NP+1),IP(MP),IN(MP,2)
C
      IF(MP.LT.NP) THEN
        WRITE(6,*) 'MP should be, MP.GE.NP'
        STOP
        ENDIF
C
      N1=N+1
      DET=1.0
      DO 5 J=1,N
    5 IP(J)=0
      DO 200 I=1,N
      T=0.0
      DO 60 J=1,N
      IF(IP(J).EQ.1) GO TO 60
      DO 59 K=1,N
      IF(IP(K).EQ.1) GO TO 59
      TTX=ABS(A(J,K))
      IF(T.GT.TTX) GO TO 59
      IR=J
      IC=K
      T=TTX
   59 CONTINUE
   60 CONTINUE
      IF(I.EQ.1) TMAX=T
      IP(IC)=1
      IF(IR.EQ.IC) GO TO 7
      DET=-DET
      DO 8 L=1,N1
      T=A(IR,L)
      A(IR,L)=A(IC,L)
    8 A(IC,L)=T
    7 IN(I,1)=IR
      IN(I,2)=IC
      PI=A(IC,IC)
      DET=DET*PI
      IF(ABS(PI)/TMAX.LT.1.E-6) GO TO 111
      A(IC,IC)=1.0
      DO 9 L=1,N1
    9 A(IC,L)=A(IC,L)/PI
      DO 100 LI=1,N
      IF(LI.EQ.IC) GO TO 100
      T=A(LI,IC)
      A(LI,IC)=0.0
      DO 11 L=1,N1
   11 A(LI,L)=A(LI,L)-A(IC,L)*T
  100 CONTINUE
  200 CONTINUE
      DO 13 I=1,N
      L=N-I+1
      IF(IN(L,1).EQ.IN(L,2))  GO TO 13
      JR=IN(L,1)
      JC=IN(L,2)
      DO 12 K=1,N
      T=A(K,JR)
      A(K,JR)=A(K,JC)
      A(K,JC)=T
   12 CONTINUE
   13 CONTINUE
      GO TO 9999
  111 WRITE(6,2000)
 2000 FORMAT(1H,'*** MATRIX A IS SINGULAR ***')
 9999 RETURN
      END      
