program new_method

implicit none

integer :: i,ir,j,k

integer, parameter :: nstrmax = 200  !##### nombre max de stream #####!
integer, parameter :: nsa = 181  !##### Nombre d'angle de diffusion #####!
integer, parameter :: nsa1 = nsa - 1

integer :: nsat1 , nn

integer :: nstr

integer :: itemp1

real(8) :: temp1,temp2,temp3,temp4

real(8) :: teta

real(8) :: dpac 

real(8), dimension(0:nsa1) :: angle,p11,p12,p22,p33,p43,p44
real(8), dimension(0:nsa1) :: t,p
real(8), dimension(0:1000) :: trac,tac,pac

real(8), dimension(0:nstrmax,0:nstrmax) :: ma,ml
real(8), dimension(0:nstrmax,0:nstrmax) :: tab
real(8), dimension(0:nstrmax,0:999) :: pl

real(8), dimension(0:nstrmax) :: mb
real(8), dimension(0:nstrmax) :: vx,vy

real(8) :: cdr = 0.01745329252  ! Conversion degree-radian.

real(8), dimension(500)  :: x1,w1,y1

real(8), dimension(0:999) :: w

real(8), dimension(0:200) :: pmom

nstr = 60  !##### nombre de stream #####!

!----------------------------------------------------------------------------------------!
! Input data

open(1,file='phase_function_01.txt')

!read(1,*) temp1,temp1,temp1
!read(1,*) temp1,temp1,temp1,temp1

do i = 0 , nsa1
read(1,*) angle(i),p11(i)!,p12(i),p22(i),p33(i),p43(i),p44(i)
enddo

close(1)

!----------------------------------------------------------------------------------------!
! Step 1 : Troncation de la fonction de phase.

p(0) = log10( 10.0d0 )  !##### Valeur de la fonction de phase à 0° #####!
dpac = 0.9d0

temp1 = 0.0
i = 1

do while( temp1 == 0.0 )

temp2 = p(0) - log10( p11(i) )

if( temp2 .ge. dpac ) then
ir = i
temp1 = 1.0
end if

i = i + 1

enddo

nsat1 = nsa - ir

j = 1

do i = ir , nsa1
t(j) = angle(i)
p(j) = log10( p11(i) )
j = j + 1
enddo

nn = nsat1

!----------------------------------------------------------------------------------------!

call gauleg(0.0d0,1.0d0,x1,w1,500)

do i = 0 , 499

trac(i)      = - x1(500-i)
trac(1001-i) =   x1(500-i)

w(i)         = w1(500-i)
w(1001-i)    = w1(500-i)

enddo

do i = 0 , 999
tac(i) = acos( trac(i) ) / cdr
enddo

call myspline(nn,t,p,1000,tac,pac)

do i = 0 , 999
pac(i) = 10 ** pac(i)
enddo

open(8,file='test.txt')

do i = 0 , 999
write(8,*) tac(i),pac(i)
enddo

close(8)

!----------------------------------------------------------------------------------------!
! Step 2 : Calcul des polynomes de Legendre.

do i = 0 , 999 

pl(0,i) = 1.0
pl(1,i) = trac(i) 

temp1 = pl(0,i)
temp2 = pl(1,i)

temp3 = trac(i)

do j = 2 , nstr 

temp4 = ( 2.0 * j - 1.0 ) * temp3 * temp2 - ( j - 1.0 ) * temp1 

pl(j,i) = ( 1.0 / j ) * temp4

temp1 = temp2
temp2 = pl(j,i)

enddo

enddo

!----------------------
! Verification.

!open(1,file='outlst_pl')

!do i = 0 , 999

!temp1 = trac(i)

!write(1,*) temp1,(pl(j,i),j=0,nstr)

!enddo

!close(1)

!----------------------------------------------------------------------------------------!
! Step 3 : Calcul des matrices A et B.

!----------------------
! Matrice A.

do i = 0 , nstr

do j = 0 , nstr

temp1 = 0.0

do k = 0 , 999

temp1 = temp1 + ( pl(i,k) * pl(j,k) ) / ( pac(k) * pac(k) )

enddo

ma(i,j) = temp1

enddo

enddo

!----------------------
! Matrice B.

do i = 0 , nstr

temp1 = 0.0

do k = 0 , 999

temp1 = temp1 + pl(i,k) / pac(k)

enddo

mb(i) = temp1 

enddo

!----------------------
! Verification.

!open(1,file='output_amatrice')
!open(2,file='output_bmatrice')

!write(1,*)

!write(1,*) 'Matrice A: hellllo'

!write(1,*)

!do i = 0 , nstr
!write(1,*) (ma(i,j),j=0,nstr)
!enddo

!write(2,*)

!write(2,*) 'Matrice B:'

!write(2,*)

!do i = 0 , nstr
!write(2,*) mb(i)
!enddo

!close(1)
!close(2)

!----------------------------------------------------------------------------------------!
! Step 4 : Factorisation de Choleski.

!----------------------
! Calcul pour j = 1.

ml(0,0) = sqrt( ma(0,0) )

do i = 1 , nstr

ml(i,0) = ma(0,i) / ml(0,0)

enddo

!----------------------
! Calcul pour j > 1.

do j = 1 , nstr

do i = 0 , nstr

    if( i .lt. j ) then

    ml(i,j) = 0.0

    else

    if( i .eq. j) then

    temp1 = 0.0

    do k = 0 , i -1
    temp1 = temp1 + ml(i,k) * ml(i,k)
    enddo

    ml(i,j) = sqrt( ma(i,i) - temp1 )

    else

    temp1 = 0.0

    do k = 0 , i - 1
    temp1 = temp1 + ml(j,k) * ml(i,k)
    enddo

    ml(i,j) = ( ma(j,i) - temp1 ) / ml(j,j)

    endif

    endif

enddo

enddo

!----------------------
! Verification.

!write(*,*)

!write(*,*) 'Verification de la decomposition de Choleski:'

!write(*,*)

!do i = 0 , nstr 

!do j = 0 , nstr

!temp1 = 0.0

!do k = 0 , nstr
!temp1 = temp1 + ml(i,k) * ml(j,k)
!enddo

!tab(i,j) = temp1

!enddo

!enddo

!do i = 0 , nstr
!write(*,*) (tab(i,j),j=0,nstr)
!enddo

!----------------------------------------------------------------------------------------!
! Step 5 : Resolution du systeme d'equations lineaire.

!----------------------
! Resolution du 1er systeme d'equations.

vy(0) = mb(0) / ml(0,0)

do i = 1 , nstr

temp1 = 0.0

do k = 0 , i - 1
temp1 = temp1 + ml(i,k) * vy(k)
enddo

vy(i) = ( mb(i) - temp1 ) / ml(i,i)

enddo

!----------------------
! Resolution du second systeme d'equations.

vx(nstr) = vy(nstr) / ml(nstr,nstr)

do i = nstr - 1 , 0 , -1

temp1 = 0.0

do k = nstr , i + 1 , -1
temp1 = temp1 + ml(k,i) * vx(k)
enddo

vx(i) = ( vy(i) - temp1 ) / ml(i,i)

enddo

!----------------------
! Ecriture de la solution.

open(1,file='res_new_method1.txt')

do i = 0 , nstr
vx(i) = vx(i) / ( 2.0d0 * i + 1.0d0 )
write(1,*) vx(i) / vx(0)
enddo

close(1)

!----------------------------------------------------------------------------------------!
! Reconstruction of the phase function from Legendre polynomials 






!----------------------------------------------------------------------------------------!

contains

subroutine gauleg(x1,x2,x,w,n)

integer :: n  
integer :: i,j,m

real(8) :: x1,x2  

real(8), dimension(n) :: x,w  

double precision, parameter :: eps = 3.d-14  

double precision :: p1,p2,p3,pp,xl,xm,z,z1

m =  ( n + 1 ) / 2
xm = 0.5d0 * ( x2 + x1 )
xl = 0.5d0 * ( x2 - x1 )

do i = 1 , m

z  = cos( 3.141592654d0 * ( i - .25d0 ) / ( n + .5d0 ) )
z1 = z + 1.0

do while( abs( z - z1 ) .gt. EPS )
        
p1 = 1.d0
p2 = 0.d0

do j = 1 , n
p3 = p2
p2 = p1
p1 = ( ( 2.d0 * j - 1.d0 ) * z * p2 - ( j - 1.d0 ) * p3 ) / j
enddo

pp = n * ( z * p1 - p2 ) / ( z * z - 1.d0 )
z1 = z
z = z1 - p1 / pp

enddo

x(i) = xm - xl * z
x(n+1-i) = xm + xl * z
w(i) = 2.d0 * xl / ( ( 1.d0 - z * z ) * pp * pp )
w(n+1-i) = w(i)

enddo

end subroutine gauleg

!----------------------------------------------------------------------------------------!

subroutine myspline( n , x , y , nn , xn , yn )

integer :: n , nn , i

real(8) :: xx , yy

real(8) :: yp1 , ypn

real(8), dimension(0:180)  :: x , y , y2
real(8), dimension(nn)  :: xn , yn 

yp1 = 1.e31
ypn = 1.e31

call spline( x , y , n , yp1 , ypn , y2 )

do i = 1 , nn
xx = xn(i)
call splint( x , y , y2 , n , xx , yy )
yn(i) = yy
enddo

end subroutine myspline

!----------------------------------------------------------------------------------------!

subroutine spline( x , y , n , yp1 , ypn , y2 )

real(8), dimension(181) :: x , y , y2 , u

real(8) :: yp1,ypn

integer :: n

real(8) :: sig , p , qn , un

if ( yp1 > 0.99e30) then
y2(1) = 0.0
u(1)  = 0.0
else
y2(1) = -0.5
u(1)  = ( 3.0 / ( x(2) - x(1) ) ) * ( ( y(2) - y(1) ) / ( x(2) - x(1) ) - yp1 )
end if
 
do i = 2 , n - 1
sig   = ( x(i) - x(i-1) ) / ( x(i+1) - x(i-1) )
p     = sig * y2(i-1) + 2.0
y2(i) = ( sig - 1.0 ) / p
u(i)  = ( 6.0 * ( ( y(i+1) - y(i) ) / ( x(i+1) - x(i) ) - ( y(i) - y(i-1) ) &
&        / ( x(i) - x(i-1) ) ) / ( x(i+1) - x(i-1) ) - sig * u(i-1) ) / p
end do    

if ( ypn > 0.99e30) then
qn = 0.0
un = 0.0
else
qn = 0.5
un = ( 3.0 / ( x(n) - x(n-1) ) ) * ( ypn - ( y(n) - y(n-1) ) / ( x(n) - x(n-1) ) )
end if

y2(n) = ( un - qn * u(n-1) ) / ( qn * y2(n-1) + 1.0 )

do k = n-1 , 1 , -1
y2(k) = y2(k) * y2(k+1) + u(k)
end do

end subroutine spline

!----------------------------------------------------------------------------------------!

subroutine splint( xa , ya , y2a , n , x , y )

real(8), dimension(181) :: xa , ya , y2a

real(8) :: x, y

integer :: klo , khi

integer :: n

real(8) :: h , b , a

klo = 1
khi = n

do while( khi - klo > 1.0 )
k = ( khi + klo ) / 2.0
if( xa( k ) > x) then
khi = k
else
klo = k
end if
end do
	
h = xa( khi ) - xa ( klo )
if( h == 0.0 ) then
write(6,*) 'probleme: bad xa input.'
else
a = ( xa( khi ) - x ) / h
b = ( x - xa( klo ) ) / h
y = a * ya( klo ) + b * ya( khi ) + &
&( (a * a * a - a ) * y2a( klo ) + ( b * b * b - b ) * y2a( khi ) ) * ( h * h ) / 6.0
end if

end subroutine splint

!----------------------------------------------------------------------------------------!

subroutine calmom(x,y,w,nstr,pmom)

real(8), dimension(0:999) :: x , y , w
real(8), dimension(0:200) :: pmom

integer :: nstr, i, j, k

real(8), dimension(0:200) :: pl
real(8), dimension(0:200,0:999) :: apl

do i = 0 , 999

call fleg(x(i),pl,101)

do j = 0 , nstr
apl(j,i) = pl(j+1)
enddo

enddo

!----------------------

do j = 0 , nstr

pmom(j) = 0.

do i = 0 , 999
pmom(j) = pmom(j) + w(i) * y(i) * apl(j,i) * 0.5	
enddo

enddo

!----------------------

if ( abs( pmom(0) - 1. ) .gt. 0.05 ) then
do i = 0 , 999
y(i) = y(i) / pmom(0)
enddo
endif

do j = 0 , nstr
pmom(j) = pmom(j) / pmom(0)
enddo
pmom(0) = 1.

end subroutine calmom

!----------------------------------------------------------------------------------------!

subroutine fleg(x,pl,nl)

integer :: nl

real(8) :: x,pl(nl)

integer :: j

real(8) :: d,f1,f2,twox

pl(1) = 1.
pl(2) = x

if( nl .gt. 2 ) then

twox = 2. * x
f2 = x
d = 1.

do j = 3 , nl
f1 = d
f2 = f2 + twox
d = d + 1.
pl(j) = ( f2 * pl(j-1) - f1 * pl(j-2) ) / d
enddo

endif

end subroutine fleg

!----------------------------------------------------------------------------------------!

end program new_method

About / Au sujet de

Ben
Created:
19 octobre 2016 17:02:04

Updated:
19 octobre 2016 17:02:04

License / Licence

MIT License