Uebung 7 (#4)

* Ue 7

* ue update

* commit to branch

* c

* final commit

Uebung 7/3a).gp

@@ -0,0 +1,36 @@
+set term pngcairo size 2000,1080
+set output '3b)phase_space.png'
+
+
+set title 'phasespace trajectories --  free pendulum'
+set xlabel 'theta [rad]'
+set ylabel 'theta dot [rad/s]'
+
+set grid
+plot 'ruku0.01.dat' u 2:3 every 2 w lines title "0.01pi", 'ruku0.1.dat' u 2:3 every 2  w lines title "0.1pi" , 'ruku0.5.dat' u 2:3 every 2  w lines title "0.5pi", 'ruku0.99.dat' u 2:3 every 2  w lines title "0.99pi"
+
+unset out
+
+unset term
+
+set term pngcairo #size 2000,1080
+set output '3b)fourier.png'
+
+set logscale y
+set format y "10^{%L}"
+set xrange [:5]
+
+set title 'Fourier Transform -- y-axis rescaled'
+set xlabel 'omega'
+set ylabel 'spectrum'
+
+set grid
+
+plot 'fourier0.01.dat' u 1:2 w lines title "0.01pi"  ,'fourier0.1.dat' u 1:(0.01*$2) w lines title "0.01pi"  ,'fourier0.5.dat' u 1:((1./50.)**2*$2) w lines title "0.5pi" ,'fourier0.99.dat' u 1:((1./99.)**2*$2) w lines title "0.99pi"
+
+#plot 'fourier0.01.dat' u 1:2 w lines  ,'fourier0.1.dat' u 1:2 w lines  ,'fourier0.5.dat' u 1:2 w lines  ,'fourier0.99.dat' u 1:2 w lines
+
+
+unset out
+
+unset term

Uebung 7/3b).gp

@@ -0,0 +1,61 @@
+set term pngcairo size 2000,1080
+set output '3b)phase_space1.png'
+#set title 'phasespace trajectories --  free pendulum'
+set xlabel 'theta [rad]'
+set ylabel 'theta dot [rad/s]'
+set grid
+plot 'ruku_eta0.01.dat' u 2:3 every 2 w lines title "eta = 0.01 1/s"
+unset out
+unset term
+
+set term pngcairo size 2000,1080
+set output '3b)phase_space2.png'
+#set title 'phasespace trajectories --  free pendulum'
+set xlabel 'theta [rad]'
+set ylabel 'theta dot [rad/s]'
+set grid
+plot 'ruku_eta0.3.dat' u 2:3 every 2 w lines title "eta = 0.3 1/s"
+unset out
+unset term
+
+set term pngcairo size 2000,1080
+set output '3b)phase_space3.png'
+#set title 'phasespace trajectories --  free pendulum'
+set xlabel 'theta [rad]'
+set ylabel 'theta dot [rad/s]'
+set grid
+plot 'ruku_eta0.5.dat' u 2:3 every 2 w lines title "eta = 0.5 1/s"
+unset out
+unset term
+
+set term pngcairo size 2000,1080
+set output '3b)phase_space4.png'
+#set title 'phasespace trajectories --  free pendulum'
+set xlabel 'theta [rad]'
+set ylabel 'theta dot [rad/s]'
+set grid
+plot 'ruku_eta0.9.dat' u 2:3 every 2 w lines title "eta = 0.9 1/s"
+unset out
+unset term
+
+
+
+
+set term pngcairo #size 2000,1080
+set output '3b)fourier.png'
+
+set logscale y
+set format y "10^{%L}"
+set xrange [:5]
+
+set title 'Fourier Transform -- y-axis rescaled'
+set xlabel 'omega'
+set ylabel 'spectrum'
+
+set grid
+
+plot 'fourier_eta0.01.dat' u 1:2 w lines title "eta = 0.01 1/s"  ,'fourier_eta0.3.dat' u 1:2 w lines title "eta = 0.3 1/s"  ,'fourier_eta0.5.dat' u 1:2 w lines title "eta = 0.5 1/s" ,'fourier_eta0.9.dat' u 1:2 w lines title "eta = 0.9 1/s"
+
+unset out
+
+unset term

Uebung 7/3c).gp

@@ -0,0 +1,38 @@
+set term pngcairo size 2000,1080
+set output '3c)phase_space.png'
+
+
+set title 'phasespace trajectories --  free pendulum'
+set xlabel 'theta [rad]'
+set ylabel 'theta dot [rad/s]'
+
+set logscale x
+
+
+set grid
+p 'ruku_f1.dat' u 2:3 every 2  w lines title "f = 1 1/(s*s)" , 'ruku_f5.dat' u 2:3 every 2  w lines title "f = 5 1/(s*s)", 'ruku_f10.dat' u 2:3 every 2  w lines title "f = 10 1/(s*s)"
+unset logscale x
+
+unset out
+
+unset term
+
+
+set term pngcairo #size 2000,1080
+set output '3c)fourier.png'
+
+set logscale y
+set format y "10^{%L}"
+set xrange [:5]
+
+set title 'Fourier Transform -- y-axis rescaled'
+set xlabel 'omega'
+set ylabel 'spectrum'
+
+set grid
+
+plot 'fourier_f1.dat' u 1:2 w lines title "f = 1 1/(s*s)"       , 'fourier_f5.dat' u 1:2 w lines title "f = 5 1/(s*s)"  ,          'fourier_f10.dat' u 1:2 w lines title "f = 10 1/(s*s)"
+
+unset out
+
+unset term

Uebung 7/functions.f90

@@ -0,0 +1,111 @@
+module functions
+  use iso_fortran_env, only : real64, int64
+  use pendulum_parameter, only : eta,natural_frequency_squared, reduced_force, driver_frequency, stepsize_time
+  implicit none
+
+  private
+  public :: runge_kutta, trapezint, return_power
+
+contains
+
+  pure function runge_kutta(y_vector, timestep)
+
+    integer (int64), intent(in) :: timestep
+    real (real64), dimension(0:2) :: runge_kutta
+    real (real64), dimension(0:2), intent(in) :: y_vector
+
+    real (real64), dimension(0:2) :: z, k1, k2, k3, k4
+
+    z = y_vector
+
+    k1 = assign_k(z,timestep)
+
+    z = y_vector + k1*stepsize_time*0.5_real64
+
+    k2 = assign_k(z,timestep)
+
+    z = y_vector + k2*stepsize_time*0.5_real64
+
+    k3 = assign_k(z, timestep)
+
+    z = y_vector + k3*stepsize_time
+
+    k4 = assign_k(z,timestep)
+
+    runge_kutta = y_vector + stepsize_time*(k1+2*k2+2*k3+k4)/6._real64
+
+  end function runge_kutta
+
+  pure function assign_k(z,timestep)
+
+    integer (int64), intent(in) :: timestep
+    real (real64), dimension(0:2), intent(in) :: z
+    real (real64), dimension(0:2) :: assign_k
+
+
+    assign_k = [1._real64, z(2),&
+         & -2._real64*eta*z(2) - &
+         & natural_frequency_squared*sin(z(1)) + &
+         & reduced_force*sin(timestep*stepsize_time*driver_frequency)]
+
+  end function assign_k
+
+
+  pure real (real64) function trapezint(array,ndim,startintervall,endintervall,ersterindex,letzterindex) result(res)
+    integer (int64) :: i
+    integer (int64), intent(in) :: ndim
+    integer (int64), optional, intent(in) :: ersterindex,letzterindex
+    real (real64), dimension(ndim), intent(in) :: array
+    real (real64), intent(in) :: startintervall,endintervall
+    integer (int64) :: startwert,endwert
+    real (real64) :: hilf, schrittweite
+
+    if ((present(ersterindex) .AND. present(letzterindex))) then
+       startwert = ersterindex
+       endwert   = letzterindex
+    else if (present(ersterindex)) then
+       startwert = ersterindex
+       endwert = ndim
+    else if (present(letzterindex)) then
+       startwert = 1_int64
+       endwert = letzterindex
+    else
+       startwert = 1_int64
+       endwert   = ndim
+    endif
+
+
+    schrittweite = (endintervall-startintervall)/ndim
+
+    hilf = 0.
+
+    do i = startwert, endwert, 1
+       if(i == startwert .OR. i == endwert) then
+          hilf = hilf + array(i)*0.5_real64
+       else
+          hilf = hilf + array(i)
+       end if
+    end do
+
+    res = hilf*schrittweite
+    !
+  end function trapezint
+
+
+
+  pure integer (int64) function return_power(arg)
+    integer (int64), intent(in) :: arg
+
+    !POPCNT(I) returns the number of bits set (’1’ bits) in the binary representation of I.
+
+    if(arg <= 0 .or. popcnt(arg) > 1) then
+       return_power = -1
+       return
+    end if
+    !TRAILZ returns the number of trailing zero bits of an integer. Fortran 2008 and later
+
+    return_power = trailz(arg)
+
+  end function return_power
+
+end module functions

Uebung 7/pendulum_parameter.f90

@@ -0,0 +1,44 @@
+module pendulum_parameter
+  !-----------------------------------------------------------
+  !  Module pendpar contains the parameters of the driven
+  !  pendulum
+  !-----------------------------------------------------------
+  use iso_fortran_env, only : real64, int64
+  implicit none
+
+  real (real64), parameter :: PI = 4 * atan (1.0_real64)
+  real (real64), parameter ::  gravitational_acceleration = 9.80665_real64
+  real (real64), parameter :: hh=6.626076e-34_real64 ! Js ! h
+
+
+  real (real64), parameter :: length_of_pendulum = 2.451662499_real64   !meter
+  real (real64), parameter :: mass_of_pendulum = 0.01_real64            !kg
+  real (real64), parameter :: eta = 0._real64                          !1/s
+  real (real64), parameter :: reduced_force = 0._real64                !1/s**2
+  real (real64), parameter :: driver_frequency = 0._real64!2.1_real64             !1/s
+
+  ! derived quantities
+
+  real (real64), parameter :: natural_frequency_squared = &
+       &gravitational_acceleration/length_of_pendulum
+  real (real64), parameter :: natural_frequency = &
+       &sqrt(gravitational_acceleration/length_of_pendulum)
+
+  real (real64), parameter :: stepsize_time = PI/(natural_frequency*30._real64)
+  real (real64), parameter :: h_bar_times_omega0 = &
+       & hh*natural_frequency/(2._real64*PI)
+
+
+
+  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 pendulum_parameter

Uebung 7/ue7bsp1.f90

@@ -0,0 +1,70 @@
+program ue7
+  use iso_fortran_env, only :  real64, real32, int64, output_unit
+  use functions, only :  runge_kutta, trapezint, return_power
+  use pendulum_parameter, only : stepsize_time,mass_of_pendulum, h_bar_times_omega0, length_of_pendulum
+
+
+  implicit none
+
+  integer (int64) ::  i,l, ioerr
+  real (real64), parameter :: PI = 4 * atan (1.0_real64)
+  integer (int64), parameter :: number_of_timesteps = 2._real64**14
+
+  real (real64), dimension(0:2) :: y
+
+  real(real32), dimension(number_of_timesteps) :: theta_values
+  real(real32), dimension(number_of_timesteps) :: theta_values_fourier
+  real(real64), dimension(number_of_timesteps) :: theta_values_real64
+  real(real64), dimension(number_of_timesteps) :: theta_values_dot_real64
+
+  !initial values
+
+  real (real64), parameter :: initial_theta = PI*0.0_real64
+  real (real64), parameter :: initial_theta_dot = sqrt(h_bar_times_omega0/(mass_of_pendulum*length_of_pendulum**2)) * length_of_pendulum  !0._real64
+
+  real (real64) ::  sampling_parameter
+
+
+  theta_values = 0
+  theta_values_fourier = 0
+  l=0
+  sampling_parameter = 2._real64*PI/(number_of_timesteps*stepsize_time)
+
+  y = [0._real64, initial_theta, initial_theta_dot]
+
+  open(21, file="ruku.dat", status='replace')
+
+
+  do i = 1, number_of_timesteps, 1
+     y = runge_kutta(y, timestep = i)
+     write(21,*)  y(0), y(1), y(2)
+     if (mod(i,number_of_timesteps/10)==0) then
+        l=l+10
+        write(output_unit,*) l, "%"
+     end if
+     theta_values(i) = real(y(1), real32)
+     theta_values_real64(i) = y(1)
+     theta_values_dot_real64 = y(2)
+  end do
+
+  close(unit=21, iostat=ioerr, status="keep")
+  if ( ioerr /= 0 ) stop "Error closing file unit 3"
+
+  call RPA(return_power(number_of_timesteps), theta_values, 1, theta_values_fourier, 1)
+  !
+  open(unit=22, file="6.7s.dat")
+  if ( ioerr /= 0 ) stop "Error opening file "
+
+  do i = 1, number_of_timesteps/2, 1
+     write(22, *) sampling_parameter*i, theta_values_fourier(i)**2 + theta_values_fourier(i+number_of_timesteps/2)**2
+  end do
+
+  close(unit=22, status="keep")
+
+  write(*,*) stepsize_time
+
+  write(*,*)  maxval(theta_values_real64), maxval(theta_values_dot_real64), &
+       &maxval(theta_values_real64)*maxval(theta_values_dot_real64)*Pi*mass_of_pendulum
+
+
+end program ue7