program main

implicit none 

integer, parameter :: pr = selected_real_kind(15,3)

integer :: i,j,k

real(pr) :: a

integer :: n

integer :: ierr, nb_lines

real(pr), dimension(:), allocatable :: angle,p11,p12,p22,p33,p43,p44

real(pr), dimension(:,:), allocatable :: pl

real(pr) :: pl_n2, pl_n1

real(pr) :: pi = dacos(-1.0d0)

real(pr) :: cdr

real(pr) :: x

real(pr), dimension(:,:), allocatable :: MA, ML
real(pr), dimension(:), allocatable :: MB, MY, MX, p11_lp

cdr = pi / 180.0

!----------------------------------------------------------------------------------------!
! user input 

n = 120 ! Number of Legendre Polynomial Coefficients 

open(1,file='hexagomal_pristine_crystal_GOM.txt', &
    & form="formatted",iostat=ierr,status="old")
	do i = 1, 18
		read(1,*)
	end do
	if (ierr/=0) then
		write(*,*) "Impossible to open the file: modis_associated_files.txt"
		stop
	endif
	nb_lines = 0
	do while (ierr==0)
		read(1,*,iostat=ierr)
		if (ierr==0) nb_lines = nb_lines + 1
	enddo
close(1)

!write(6,*) 'nb_lines', nb_lines

allocate( angle( nb_lines ) )
allocate( p11( nb_lines ) )
allocate( p12( nb_lines ) )
allocate( p22( nb_lines ) )
allocate( p33( nb_lines ) )
allocate( p43( nb_lines ) )
allocate( p44( nb_lines ) )

open(1,file='hexagomal_pristine_crystal_GOM.txt')
	do i = 1, 18
		read(1,*)
	end do
	do i = 1 , nb_lines
		read(1,*) angle(i),p11(i),p12(i),p22(i),p33(i),p43(i),p44(i)
		!write(6,*) angle(i)
	enddo
close(1)

!----------------------------------------------------------------------------------------!
! Legendre Polynomials 

allocate( pl( n, nb_lines ) )

do i = 1, nb_lines
	x = dcos( angle(i) * cdr )
    pl(1,i) = 1.0
    pl(2,i) = x
    pl_n2 = pl(1,i)
    pl_n1 = pl(2,i)
    do j = 3 , n
        pl(j,i) = ( 1.0 / j ) * ( ( 2.0 * j - 1.0 ) * x * pl_n1 - ( j - 1.0 ) * pl_n2 )
        pl_n2 = pl_n1
        pl_n1 = pl(j,i)
    end do
end do

!----------------------------------------------------------------------------------------!
! A matrix

allocate( MA(n,n) )

do i = 1 , n
	do j = 1 , n
	a = 0.0
		do k = 1, nb_lines
			a = a + ( pl(i,k) * pl(j,k) ) / ( p11(k)**2 )
		enddo
	MA(i,j) = a
	enddo
enddo

!----------------------------------------------------------------------------------------!
! B matrix

allocate( MB(n) )

do i = 1 , n
	a = 0.0
	do k = 1, nb_lines
		a = a + pl(i,k) / p11(k)
	enddo
	MB(i) = a 
enddo

!----------------------------------------------------------------------------------------!
! Choleski Factorisation 

allocate( ML(n,n) )

ML(1,1) = sqrt( MA(1,1) )
do i = 2 , n
	ML(i,1) = MA(1,i) / ML(1,1)
enddo

do j = 2 , n
	do i = 1 , n
		if( i .lt. j ) then
			ML(i,j) = 0.0
		else
			if( i .eq. j) then
				a = 0.0
				do k = 1 , i -1
					a = a + ML(i,k)**2
				enddo
				ML(i,j) = sqrt( MA(i,i) - a )
			else
				a = 0.0
				do k = 1 , i - 1
					a = a + ML(j,k) * ML(i,k)
				enddo
				ML(i,j) = ( MA(j,i) - a ) / ML(j,j)
			endif
		endif
	enddo
enddo

!----------------------------------------------------------------------------------------!
! Solve linear system 1

allocate( MY(n) )

MY(1) = MB(1) / ML(1,1)
do i = 2 , n
	a = 0.0
	do k = 1 , i-1
		a = a + ML(i,k) * MY(k)
	enddo
	MY(i) = ( MB(i) - a ) / ML(i,i)
enddo

!----------------------------------------------------------------------------------------!
! Solve linear system 2

allocate( MX(n) )

MX(n) = MY(n) / ML(n,n)
do i = n-1, 1, -1
	a = 0.0
	do k = n, i+1, -1
		a = a + ML(k,i) * MX(k)
	enddo
	MX(i) = ( MY(i) - a ) / ML(i,i)
enddo

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

allocate( p11_lp(nb_lines) )

do i = 1, nb_lines
	a = 0
    do j = 1 , n
    	a = a + MX(j) * pl(j,i) 
    end do
	p11_lp(i) = a
end do

!----------------------------------------------------------------------------------------!
! Save results

open(1,file='phase_function_legendre_polynomial_expansion_01.txt')
do i = 1 , n
	MX(i) = MX(i) / ( 2.0d0 * i + 1.0d0 )
	write(1,*) MX(i) / MX(1)
	!write(6,*) MX(i) !/ MX(1)
enddo
close(1)

open(1,file='phase_function_legendre_polynomial_expansion_02.txt')
do i = 1 , nb_lines
	write(1,*) angle(i), p11_lp(i), p11(i)
enddo
close(1)

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

end program main

Infos

Main author / Auteur principal:
Ben

Created / Créé: 2016-10-19

Edited / Modifié: 2016-10-19

License / Licence: MIT License

Abstract / Résumé