program Fortran_Figure6_6_a

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: N=3000
 integer, parameter :: Npom=500
 integer, parameter :: Nsaves=20
 real(kind=long), parameter :: pi=3.14159265358979_long
 real(kind=long), dimension(3,Nsaves) :: Ksave 
 real(kind=long), dimension(Nsaves) :: Psave
 real(kind=long), dimension(N,N) :: A,K
 real(kind=long), dimension(N) ::  rhs,r,g
 real(kind=long) ::  gam
 real(kind=long) ::  b,h1,h2,plambda,pom,pom2,pom3,pom4,knum,ri,Ki,intKg
 integer :: ii,jj,kk,ns,no
 

 do ns=1,Nsaves
    Psave(ns)=real(ns)/20_long
    do ii=1,3
       Ksave(ii,ns)=0.0_long
    end do
 end do

 do no=1,3
 
    if (no.EQ.1) then
       gam=0.5_long
    else
       if (no.EQ.2) then
          gam=1.0_long
       else
          gam=2.0_long
       end if
    end if
     
    write(*,*) 'Computation for gamma =',gam

    do ns=1,Nsaves

       plambda=Psave(ns)
       b=gam*200.0_long
       h1=(b-1.0_long)/real(N-Npom)
       h2=1.0_long/real(Npom)

       do ii=1,Npom
          r(ii)=h2/2.0_long+h2*real(ii-1)
       end do     
       do ii=Npom+1,N
          r(ii)=1.0_long+h1/2.0_long+h1*real(ii-1-Npom)
       end do     
      
       pom=sqrt(2.0_long*pi)*gam
       pom2=2.0_long*gam*gam
       pom3=sqrt(2.0_long)*gam
     
       do ii=1,N
          do jj=1,Npom
             K(ii,jj)=(r(jj)/(r(ii)*pom))*(exp(-(r(ii)-r(jj))*(r(ii)-r(jj))/(pom2))-exp(-(r(ii)+r(jj))*(r(ii)+r(jj))/(pom2)))
             A(ii,jj)=-K(ii,jj)*h2*(1.0_long-plambda)
          end do
          do jj=Npom+1,N
             K(ii,jj)=(r(jj)/(r(ii)*pom))*(exp(-(r(ii)-r(jj))*(r(ii)-r(jj))/(pom2))-exp(-(r(ii)+r(jj))*(r(ii)+r(jj))/(pom2)))
             A(ii,jj)=-K(ii,jj)*h1
          end do
          rhs(ii)=(gam*gam*K(ii,N)/b)+1_long-0.5_long*erf((b-r(ii))/(pom3))-0.5_long*erf((r(ii)+b)/(pom3))
          A(ii,ii)=1.0_long+A(ii,ii)
       end do   
      
       do ii=1,N-1
          do jj=ii+1,N
             pom4=A(jj,ii)/A(ii,ii)
             do kk=ii,N
	        A(jj,kk)=A(jj,kk)-A(ii,kk)*pom4
	     end do
	     rhs(jj)=rhs(jj)-rhs(ii)*pom4
          end do	     
       end do
        
       do ii=1,N
          pom4=A(ii,ii)
          do jj=ii,N
             A(ii,jj)=A(ii,jj)/pom4
          end do
          rhs(ii)=rhs(ii)/pom4
       end do
       
       do ii=2,N
          do jj=1,ii-1
             pom4=A(jj,ii)
	     do kk=ii,N
	        A(jj,kk)=A(jj,kk)-pom4*A(ii,kk)   	  
             end do
	     rhs(jj)=rhs(jj)-pom4*rhs(ii)
          end do
       end do
       
       do ii=1,N
          g(ii)=rhs(ii)
       end do   
          	  
       knum=0.0_long
 
       do ii=1,Npom
          knum=knum+h2*4.0_long*pi*r(ii)*r(ii)*g(ii)
       end do
 
       Ksave(no,ns)=knum*plambda
    
       write(*,*) 'Computed ns =',ns,'out of 20'

    end do
         
    open(unit=1, file='data_Figure6_6_a.dat', status='replace')

    100 format(e12.5,x,e12.5,x,e12.5,x,e12.5,x,e12.5,x,e12.5)

    write(1,100) 0.0_long,0.0_long,0.0_long,0.0_long
    do ii=1,Nsaves
       write(1,100) Ksave(1,ii),Ksave(2,ii),Ksave(3,ii),Psave(ii)
    end do

    close(unit=1)

 end do
     
end program Fortran_Figure6_6_a