Skip to content

Commit

Permalink
Uebung 8 (#5)
Browse files Browse the repository at this point in the history 
* abgabe
  • Loading branch information
@saasaa
saasaa authored May 22, 2017
1 parent 44614fb commit ff2fe86
Show file tree
Hide file tree
Showing 6 changed files with 400 additions and 0 deletions.
24 changes: 24 additions & 0 deletions Uebung 8/argand.gp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
set term pdfcairo #enhanced
set out 'argand.pdf'

#set sample 10000
unset key
set grid

set title "V_1=0, l=2"


set xlabel "Re[S]"
set ylabel "Im[S]"

plot "argand_l=2_v1=0.dat" u 1:2 w lines

set title "V_1=18 MeV, l=2"

plot "argand_l=2_v1=18.dat" u 1:2 w lines


unset out
unset term

#plot 'data.dat' index i matrix with image using 1:3 every 1:999:1:2
35 changes: 35 additions & 0 deletions Uebung 8/delta.gp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
set term pdfcairo #enhanced
set out 'delta.pdf'

#set sample 10000
unset key
set grid

set title "V_1=0, l=2"

set logscale x

set xlabel "absoute Streuphase [rad]"
set ylabel "Energie [MeV]"

plot "delta_l=2_v1=0.dat" u 1:2 w lines

set title "V_1=18 MeV, l=2"


plot "delta_l=2_v1=18.dat" u 1:2 w lines


set title "V_1=0 MeV, l=2, Streuphase"

plot "streu_delta_l=2_v1=0.dat" u 1:2 w lines

set title "V_1=18 MeV, l=2, Streuphase"

plot "streu_delta_l=2_v1=18.dat" u 1:2 w lines


unset out
unset term

#plot 'data.dat' index i matrix with image using 1:3 every 1:999:1:2
157 changes: 157 additions & 0 deletions Uebung 8/functions_ue8.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,157 @@
module functions_ue8
use iso_fortran_env, only : real64, int64
use parameter
implicit none

! private
! public :: runge_kutta, trapezint, return_power

contains

pure real (real64) function Vpot(x)
real (real64), intent(in) :: x

Vpot=-V0/(1._real64+exp((x-R0)/a0)) + V1*exp(-(x-R1)**2/a1)

end function Vpot

pure real (real64) function w(x,l,energy)
real (real64), intent(in) :: x, energy
integer (int64), intent(in) :: l

w = mhb2*energy-mhb2*Vpot(x) &
&-l*(l+1)/(x**2)

end function w

function numerov(energy, l, direction_forward)
integer (int64) :: i
logical, optional, intent(in) :: direction_forward
integer (int64), intent(in) :: l
real (real64), dimension(0:Nrmax) :: numerov
real (real64), intent(in) :: energy

real (real64), dimension(0:Nrmax) :: u
real (real64) :: q_plus_one, q, q_minus_one

u=0._real64

if (present(direction_forward) .and. direction_forward == .true. .or. (.not. present(direction_forward)) ) then

u(0) = 0._real64
u(1) = 0.01_real64

q_minus_one = 0._real64
q = (1._real64+&
&stepsize_r*stepsize_r*&
&w(stepsize_r*1._real64,l,&
&energy)/12._real64)*u(1)

do i = 2, Nrmax, 1
q_plus_one = 12._real64 * u(i-1) - 10._real64*q - q_minus_one
u(i) = q_plus_one/(1._real64+&
&stepsize_r*stepsize_r*w(stepsize_r*i,l,&
&energy)/12._real64)

q_minus_one = q
q = q_plus_one

end do

else if (present(direction_forward) .and. direction_forward == .false. ) then

u(Nrmax) = exp(-sqrt(mhb2*energy)*Nrmax*stepsize_r)
u(Nrmax-1) = exp(-sqrt(mhb2*energy)*(Nrmax-1)*stepsize_r)

q_plus_one = (1._real64+&
&stepsize_r*stepsize_r*&
&w(stepsize_r*Nrmax,l,&
&energy)/12._real64)*u(Nrmax)

q = (1._real64+&
&stepsize_r*stepsize_r*&
&w(stepsize_r*(Nrmax-1),l,&
&energy)/12._real64)*u(Nrmax-1)

write(*,*) u(Nrmax), u(Nrmax-1), q_plus_one, q


do i = 2, Nrmax-1
q_minus_one = 12._real64 * u(Nrmax-i+1) - 10._real64*q - q_plus_one
u(i) = q_minus_one/(1._real64+&
&stepsize_r*stepsize_r*w(stepsize_r*(Nrmax-i),l,&
&energy)/12._real64)

q_plus_one = q
q = q_minus_one

end do



end if

numerov = u

end function numerov


pure subroutine bessel(x,l,bes,neu)
!
! evaluates the Bessel and Neumann functions in
! Riccati form for given l at the argument x
! j0=sin(x),n0=cos(x)
! j1=sin(x)/x-cos(x)
! n1=cos(x)/x+sin(x)
! j2=sin(x)*(3/x**2-1)-3*cos(x)/x
! n2=cos(x)*(3/x**2-1)+3*sin(x)/x
!
!
integer (int64), intent(in) :: l
real (real64), intent(in) :: x
real (real64), intent(out) :: bes,neu
!
bes=0._real64
neu=0._real64
if(l.eq.0) then
bes=sin(x)
neu=cos(x)
end if
!
if(l.eq.1) then
bes=sin(x)/x-cos(x)
neu=cos(x)/x+sin(x)
end if
!
if(l.eq.2) then
bes=sin(x)*(3/x**2-1._real64)-3._real64*cos(x)/x
neu=cos(x)*(3/x**2-1._real64)+3._real64*sin(x)/x
end if

return
end subroutine bessel


end module functions_ue8

! u(Nrmax) = exp(-sqrt(mhb2*energy)*Nrmax*stepsize_r)
! u(Nrmax-1) = exp(-sqrt(mhb2*energy)*(Nrmax-1)*stepsize_r)
!
! q_plus_one = (1._real64+&
! &stepsize_r*stepsize_r*&
! &w(stepsize_r*Nrmax,l,&
! &energy)/12._real64)*u(Nrmax)
!
! q = (1._real64+&
! &stepsize_r*stepsize_r*&
! &w(stepsize_r*(Nrmax-1),l,&
! &energy)/12._real64)*u(Nrmax-1)
!
! do i = Nrmax-2, 0, -1
! q_minus_one = 12._real64 * u(i+1) - 10._real64*q - q_plus_one
! u(i) = q_minus_one/(1._real64+&
! &stepsize_r*stepsize_r*w(stepsize_r*i,l,&
! &energy)/12._real64)
!
! q_plus_one = q
! q = q_minus_one
66 changes: 66 additions & 0 deletions Uebung 8/parameter.f90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@
module parameter
!-----------------------------------------------------------
! Module pendpar contains the parameters of the driven
! pendulum
!-----------------------------------------------------------
use iso_fortran_env, only : real64, int64
implicit none

integer (int64), parameter :: Nrmax = 2000
integer (int64), parameter :: NEmax = 40000

real (real64), parameter :: PI = 4 * atan (1.0_real64)
real (real64), parameter :: hbar=6.582122e-22_real64 ! MeVs
real (real64), parameter :: hbaratto=6.582122e-4_real64 ! MeVas



real (real64), parameter :: mass_of_neutron=939.56563_real64 ! MeV
real (real64), parameter :: atomic_mass_unit=931.49432_real64 ! MeV
real (real64), parameter :: hbar_times_c=197.327053_real64 ! MeVfm

integer (int64), parameter :: A = 40 ! mass number of target
real (real64), parameter :: mA = 39.962591_real64 ! rel. atomic mass of target nucleus
real (real64), parameter :: Anorm = 1._real64 ! start normalization of wave function
real (real64), parameter :: energy_max = 1000._real64 ! maximal energy in MeV
real (real64), parameter :: radius_max = 30._real64
integer (int64), parameter :: number_of_meshpoints = 2000

!Potential parameters
real (real64), parameter :: V0 = 84._real64 ! MeV
real (real64), parameter :: rr = 1.20_real64 !fm
real (real64), parameter :: a0 = 0.6_real64 !fm

! parameters of additiona Gauss potential
!real (real64), parameter :: V1 = 18._real64 !MeV
real (real64), parameter :: V1 = 0._real64 !MeV
real (real64), parameter :: R1 = 6.00_real64 !fm
real (real64), parameter :: a1 = 4.0_real64 !fm**2

!argand
complex (real64), parameter :: ci=(0,1)
complex (real64), dimension(NEmax) :: S
real (real64), parameter :: alpha = 0.02_real64

! derived quantities
real (real64), parameter :: reduced_mass = &
&mass_of_neutron*mA*atomic_mass_unit/(mA*atomic_mass_unit+mass_of_neutron) !MeV
real (real64), parameter :: mhb2 = 2._real64*reduced_mass/(hbar_times_c**2)

real (real64), parameter :: energy_step = energy_max/NEmax
real (real64), parameter :: stepsize_r = radius_max/Nrmax
real (real64), parameter :: R0 = rr*A**(1._real64/3._real64)


! private
! public :: eta
! public :: reduced_force
! public :: driver_frequency
! public :: natural_frequency_squared
! public :: stepsize_time
! public :: mass_of_pendulum
! public :: natural_frequency
! public :: h_bar_times_omega0
! public :: length_of_pendulum

end module parameter
24 changes: 24 additions & 0 deletions Uebung 8/timedelay.gp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
set term pdfcairo #enhanced
set out 'time_delay.pdf'

#set sample 10000
unset key
set grid

set title "timedelay: V_1=18, l=2"


set xlabel "Energy [MeV]"
set ylabel "time [as]"

set logscale xy

plot "time_delay_v1=18.dat" u 1:2 w lines




unset out
unset term

#plot 'data.dat' index i matrix with image using 1:3 every 1:999:1:2
Loading

0 comments on commit ff2fe86

Please sign in to comment.