program rectangular_aperture

implicit none

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

integer :: i,j

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

real(pr) :: wl ! WaveLength (micrometer) !
real(pr) :: wk ! WaveNumber (micrometer) !
real(pr) :: lx , ly ! aperture sides !

real(pr) :: D ! Distance from the screen to the aperture (micrometer) !
real(pr) :: x_min , x_max ! screen x_min , x_max !
real(pr) :: y_min , y_max ! screen y_min , y_max !
real(pr) :: resol_x , resol_y ! x and y resolution !

integer :: nb_x_point , nb_y_point

real(pr) :: tmp_01 , tmp_02 , tmp_03 , tmp_04
real(pr) :: x , y
real(pr) :: r
real(pr) :: Intensity

integer :: output ! if 1 -> 3D !
! if 2 -> 2D !

!----------------------------------------------------------------------------------------!
! User Inputs

wl = 0.55d0
lx = 10.0d0
ly = 10.0d0
D = 1000000.0d0
x_min = - 300000.0d0
x_max =	  300000.0d0
y_min = - 300000.0d0
y_max =   300000.0d0
resol_x = 1000.0d0
resol_y = 1000.0d0
output = 1

!----------------------------------------------------------------------------------------!
! Check Resolution

wk = ( 2.0d0 * pi ) / wl

if( resol_x > ( x_max - x_min ) / 2.0d0 )then
write(6,*) "the resolution along x is too large"
end if

if( resol_y > ( y_max - y_min ) / 2.0d0 )then
write(6,*) "the resolution along y is too large"
end if

!----------------------------------------------------------------------------------------!
! Main code

open(1,file='output_data')

!----------------------------------------------------------------------------------------!
! 3d

if( output == 1 )then
nb_x_point = int( ( x_max - x_min ) / resol_x ) + 1
nb_y_point = int( ( y_max - y_min ) / resol_y ) + 1
x = x_min
do i = 1 , nb_x_point
y = y_min
do j = 1 , nb_y_point
r = dsqrt( x * x + y * y + D * D )
tmp_01 = wk * ( x / r ) * ( lx / 2.0d0 )
tmp_02 = fsc( tmp_01 )
tmp_03 = wk * ( y / r ) * ( ly / 2.0d0 )
tmp_04 = fsc( tmp_03 )
Intensity = ( tmp_02 * tmp_02 ) * ( tmp_04 * tmp_04 )
if( Intensity <= 0.01d0 )then ! To show the existence of the weak secondary maxima !
write(1,*) x , y , Intensity
else
write(1,*) x , y , 0.01d0
end if
y = y + resol_y
end do
x = x + resol_x
end do
end if

!----------------------------------------------------------------------------------------!
! 2d

if( output == 2 )then
nb_x_point = int( ( x_max - x_min ) / resol_x ) + 1
nb_y_point = int( ( y_max - y_min ) / resol_y ) + 1
x = x_min
y = 0.0d0
do i = 1 , nb_x_point
r = dsqrt( x * x + y * y + D * D )
tmp_01 = wk * ( x / r ) * ( lx / 2.0d0 )
tmp_02 = fsc( tmp_01 )
tmp_03 = wk * ( y / r ) * ( ly / 2.0d0 )
tmp_04 = fsc( tmp_03 )
Intensity = ( tmp_02 * tmp_02 ) * ( tmp_04 * tmp_04 )
write(1,*) x , y , Intensity
x = x + resol_x
end do
end if

close(1)

contains

real function fsc(x)

real(pr) :: x

if( abs(x) > 0.001d0 ) then
fsc = sin( x ) / x
else
fsc = 1.0d0 - x * x / 6.0d0
end if

return

end function fsc

end program rectangular_aperture


#### Announcement

Bienvenue sur "The Open Science Notebook Project" ! Le site est en développement, vous pouvez nous contacter par le biais du formulaire de contact si vous avez des questions ou si vous avez simplement découvert une erreur !

Date de la prochaine mise à jour le: 27 janvier 2018