program radtransf                            !Nov 29 2016
       
        
!...... This program calculates the solution for the radiative transfer
!...... equation in the Sn approximation, with n=2. The program starts 
!...... by building a Matrix for the discretized system of coupled 
!...... differential equations, obtained by the Diamond difference scheme.


!...... Version 0.0.1

        implicit none
        
        real*8 Dr,rmax
        integer npoints
        common/bkinput/ Dr,rmax,npoints


!....... Input file
        open (unit=10,file='radtrans.inp',status='old')


!....... call input        
        call readinput     

!....... Matrix generation and resolution  
        call matrix


!...... Close files

        close (unit=10)
        close (unit=17)     
        close (unit=18)   
        close (unit=19) 
        close (unit=20) 
         
        stop
        end
        
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        subroutine readinput
        implicit none     
        
!...... Input parameters

           
        real*8 Dr,rmax
        integer npoints
        common/bkinput/ Dr,rmax,npoints
        
        real*8 rzero, one, two, five
        data rzero,one,two,five/0.0d0,1.0d0,2.0d0,5.0d0/
        
        
        namelist /griddata/Dr,npoints,rmax
                
        Dr=rzero
        npoints=0
        rmax=rzero
        
!...... read griddata

        read(10,griddata)
        if (rmax.eq.rzero) rmax=Dr*npoints
        if (Dr.eq.0) Dr=rmax/npoints
        write(*,111) rmax,npoints,Dr
111     format(/5x,'rmax=',f9.3,5x,'npoints=',i8,5x,'Dr=',1pg12.5)


        return
        end

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


      subroutine matrix

!.... This program constructs a matrix A, a vector b, and then solves 
!.... the system Ax=b using the lapack package.

      
      implicit none


      real*8 Dr,rmax
      integer npoints
      common/bkinput/ Dr,rmax,npoints
      
      integer mxpts
      parameter (mxpts=10000)
      integer i,j,rows,hrows,cols,k,aux,pivot(mxpts),info

      
      real*8 A(mxpts,mxpts),b(mxpts),v(mxpts),c,x,mu1,mu2,sumaux,rsum


      real*8 one, two, four, hundred
      data one, two, four, hundred/ 1.0d0, 2.0d0, 4.0d0, 100.0d0/  


!..... Output files

!      open (unit=16,file='rhs.dat',status='unknown')
      open (unit=17,file='matrix.dat',status='unknown')
      open (unit=18,file='Psi1.dat',status='unknown')
      open (unit=19,file='Psi2.dat',status='unknown')
      open (unit=20,file='check.dat',status='unknown')

!..... Inputs
      
      rows=2*npoints
      hrows=npoints
      cols=rows !square matrix

      print*,'Rows,hrows=',rows,hrows

!..... Coefficients
!..... Initialization

        mu1=dsqrt(1.0d0/3.0d0)
        mu2=-1.0d0*dsqrt(1.0d0/3.0d0)

      do i=1,rows
        b=0.0d0
        do j=1,cols
          A(i,j)=0.0d0
       enddo
      enddo

!..... Build matrix A

      do i=1,hrows
        do j=1,cols
        if(i.eq.1 .and. i.ne.j) then
          A(i,j)=0.0d0
        else if(i.eq.1 .and. i.eq.j) then
          A(i,j)=1.0d0
        else if(i.ne.1 .and. i.eq.j) then
          A(i,j)=two*mu1
        else if(i.ne. 1 .and. i-1 .eq. j) then
          A(i,j)=-two*mu1
        else if(i.ne.1 .and. i+j.eq.rows+1) then
          A(i,j)=-hundred*Dr
        else if(i.ne.1) then
          A(i,j)=0.0d0
        endif
       enddo
      enddo

      do i=hrows+1,rows
        do j=1,cols
        if(i.eq.hrows+1 .and. i.ne.j) then
          A(i,j)=0.0d0
        else if(i.eq.hrows+1 .and. i.eq.j) then
          A(i,j)=1.0d0
        else if(i.ne.hrows+1 .and. i.eq.j) then
          A(i,j)=-two*mu2
        else if(i.ne.hrows+1 .and. i-1.eq.j) then
          A(i,j)=two*mu2
        else if(i.ne.hrows+1 .and. i+j.eq.rows+1) then
          A(i,j)=-hundred*Dr
        else if(i.ne.hrows+1) then
          A(i,j)=0.0d0
        endif
       enddo
      enddo

!        do i=1,rows
!           do j=1,cols
!                print 200,i,j,A(i,j)
!           enddo
!       enddo

!..... Vector preparation for the right hand side of Ax=b.

        do i=1,rows
                 b(i)=0.01d0*Dr
        enddo 

       b(1)=0.0d0
       b(hrows+1)=0.0d0


!	do i=1,rows
!	   print*,b(i)
!	end do


!200   format(4x,"A(",i5,",",i5,")=",1pg15.8)  !,f10.8)


!..... Print Matrix

        do i=1,rows
                do j=1,cols
                    write (17,200)i,j,A(i,j)
                enddo
        enddo

200   format(4x,"A(",i5,",",i5,")=",1pg15.8)  !,f10.8)

!      do i=1,rows
!        print*,( A(i,j), j=1,cols)
!      enddo

!..... Find the solution x=A-1*b

        call dgesv (rows, 1, A, mxpts, pivot, b, mxpts, info)

        print*,"Info=",info


	do j=1,hrows+1
           x=(j-1)*Dr
           write(18,*)x,b(j)
	end do

        
	do i=hrows,rows
           x=(i-hrows)*Dr
           write(19,*)x,b(hrows+1+rows-i)
	end do

!..... Verification for Ax=b
!..... Initialization

      do i=1,rows
        v(i)=0.0d0
      enddo

      rsum=0.0d0
      sumaux=0.0d0

!..... Verification for Ax=b

      do i=1,rows
        do j=1,cols
          sumaux=A(i,j)*b(j)
          rsum=rsum+sumaux
          v(i)=rsum
       enddo
       write(20,*)v(i)
      enddo



!      do i=1,rows
!        b(i)=0.01d0*Dr
!        write(20,*)v(i),b(i)
!      enddo 

       return    
       end

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%