mse.f

Go to the documentation of this file.

!*******************************************************************************
!
!  Note separating the Doxygen comment block here so detailed decription is
!  found in the Module not the file.
!
!*******************************************************************************
 
      MODULE mse
 
      USE stel_kinds, only: rprec, dp
      USE coordinate_utilities
      USE profiler
      USE mpi_inc
      USE signal
 
      IMPLICIT NONE
 
!*******************************************************************************
!  mse module parameters
!*******************************************************************************
      INTEGER, PARAMETER :: mse_degrees_flag = 1
      INTEGER, PARAMETER :: mse_ratio_flag   = 2
 
!*******************************************************************************
!  DERIVED-TYPE DECLARATIONS
!  1) mse base class
!
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      TYPE, EXTENDS(signal_class) :: mse_class
         REAL (rprec), DIMENSION(3) :: point
 
         REAL (rprec)               :: alpha
         REAL (rprec)               :: omega
         REAL (rprec)               :: delta
         REAL (rprec)               :: theta
 
         INTEGER                    :: flags
      CONTAINS
         PROCEDURE                  ::                                         &
     &      get_modeled_signal_last => mse_get_modeled_signal
         PROCEDURE                  ::                                         &
     &      get_type => mse_get_type
         PROCEDURE                  :: get_header => mse_get_header
         PROCEDURE, PRIVATE         ::                                         &
     &      angle_to_bt_proj => mse_angle_to_bt_proj
         final                      :: mse_destruct
      END TYPE
 
!*******************************************************************************
!  INTERFACE BLOCKS
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      INTERFACE mse_class
         MODULE PROCEDURE mse_construct_vec,                                   &
     &                    mse_construct_rad
      END INTERFACE
 
      PRIVATE :: mse_angle_to_bt_proj, mse_angle_to_horizontal
 
      CONTAINS
!*******************************************************************************
!  CONSTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION mse_construct_vec(point, view_start, view_end,                  &
     &                           beam_start, beam_end, in_degrees,             &
     &                           is_ratio)
 
      IMPLICIT NONE
 
!  Declare Arguments
      class(mse_class), POINTER  :: mse_construct_vec
      REAL (rprec), DIMENSION(3), INTENT(in) :: point
      REAL (rprec), DIMENSION(3), INTENT(in) :: view_start
      REAL (rprec), DIMENSION(3), INTENT(in) :: view_end
      REAL (rprec), DIMENSION(3), INTENT(in) :: beam_start
      REAL (rprec), DIMENSION(3), INTENT(in) :: beam_end
      LOGICAL, INTENT(in)                    :: in_degrees
      LOGICAL, INTENT(in)                    :: is_ratio
 
!  local variables
      REAL (rprec)                           :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      ALLOCATE(mse_construct_vec)
 
      mse_construct_vec%flags = 0
      IF (in_degrees) THEN
         mse_construct_vec%flags = ibset(mse_construct_vec%flags,              &
     &                                   mse_degrees_flag)
      END IF
      IF (is_ratio) THEN
         mse_construct_vec%flags = ibset(mse_construct_vec%flags,              &
     &                                   mse_ratio_flag)
      END IF
 
      mse_construct_vec%point = point
 
!  Get the angles with respect to the toroidal field component.
!  The mse alpha is defined so that 0 points in the negative torodial direction
!  and 90 points in the negative radial direction. So this needs to be PI minus
!  the angle obtained from the dotproduct.
      mse_construct_vec%alpha = pi -                                           &
     &   mse_construct_vec%angle_to_bt_proj(beam_end - beam_start)
      mse_construct_vec%omega =                                                &
     &   mse_construct_vec%angle_to_bt_proj(view_end - view_start)
 
!  Get the angles with respect to the horizontal.
!  The mse delta is defined so that 90 points in the negative z direction. So
!  this needs to be the negation of the angle obtained from the dotproduct.
      mse_construct_vec%delta =                                                &
     &   -mse_angle_to_horizontal(beam_end - beam_start)
      mse_construct_vec%theta =                                                &
     &   mse_angle_to_horizontal(view_end - view_start)
 
      CALL profiler_set_stop_time('mse_construct_vec', start_time)
 
      END FUNCTION
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION mse_construct_rad(point, alpha, omega, delta, theta,            &
     &                           in_degrees, is_ratio)
 
      IMPLICIT NONE
 
!  Declare Arguments
      class(mse_class), POINTER             :: mse_construct_rad
      REAL (rprec), DIMENSION(3), INTENT(in) :: point
      REAL (rprec), INTENT(in)               :: alpha
      REAL (rprec), INTENT(in)               :: omega
      REAL (rprec), INTENT(in)               :: delta
      REAL (rprec), INTENT(in)               :: theta
      LOGICAL, INTENT(in)                    :: in_degrees
      LOGICAL, INTENT(in)                    :: is_ratio
 
!  local variables
      REAL (rprec)                           :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      ALLOCATE(mse_construct_rad)
 
      mse_construct_rad%flags = 0
      IF (in_degrees) THEN
         mse_construct_rad%flags = ibset(mse_construct_rad%flags,              &
     &                                   mse_degrees_flag)
      END IF
      IF (is_ratio) THEN
         mse_construct_rad%flags = ibset(mse_construct_rad%flags,              &
     &                                   mse_ratio_flag)
      END IF
 
      mse_construct_rad%point = point
 
!  Get the angles with respect to the toroidal field component.
      mse_construct_rad%alpha = alpha
      mse_construct_rad%omega = omega
      mse_construct_rad%delta = delta
      mse_construct_rad%theta = theta
 
      CALL profiler_set_stop_time('mse_construct_rad', start_time)
 
      END FUNCTION
 
!*******************************************************************************
!  DESTRUCTION SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE mse_destruct(this)
 
      IMPLICIT NONE
 
!  Declare Arguments
      TYPE (mse_class), INTENT(inout) :: this
 
!  Start of executable code
      this%point = 0.0
      this%alpha = 0.0
      this%omega = 0.0
      this%delta = 0.0
      this%theta = 0.0
      this%flags = 0
 
      END SUBROUTINE
 
!*******************************************************************************
!  GETTER SUBROUTINES
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION mse_get_modeled_signal(this, a_model, sigma, last_value)
      USE model
 
      IMPLICIT NONE
 
!  Declare Arguments
      REAL (rprec), DIMENSION(4) :: mse_get_modeled_signal
      CLASS (mse_class), INTENT(inout)        :: this
      TYPE (model_class), POINTER             :: a_model
      REAL (rprec), DIMENSION(4), INTENT(out) :: sigma
      REAL (rprec), DIMENSION(4), INTENT(in)  :: last_value
 
!  local variables
      REAL (rprec), DIMENSION(3)              :: b_cyl
      REAL (rprec)                            :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      sigma = 0.0
 
      IF (btest(a_model%state_flags, model_state_vmec_flag)   .or.             &
     &    btest(a_model%state_flags, model_state_siesta_flag) .or.             &
     &    btest(a_model%state_flags, model_state_shift_flag)  .or.             &
     &    btest(a_model%state_flags, model_state_signal_flag)) THEN
 
!  Get the magnetic field vector.
         b_cyl = equilibrium_get_b_vec(a_model%equilibrium, this%point,        &
     &                                 .true.)
 
         mse_get_modeled_signal(2) =                                           &
     &         -(b_cyl(3)*cos(this%delta)*cos(this%omega + this%alpha) +       &
     &           sin(this%delta)*(b_cyl(1)*sin(this%omega) -                   &
     &                            b_cyl(2)*cos(this%omega)))
 
         mse_get_modeled_signal(3) =                                           &
     &         -b_cyl(3)*sin(this%theta)*cos(this%delta)*                      &
     &          sin(this%omega +  this%alpha)                                  &
     &       + sin(this%theta)*sin(this%delta)*                                &
     &         (b_cyl(2)*sin(this%omega) + b_cyl(1)*cos(this%omega))           &
     &       + cos(this%theta)*cos(this%delta)*                                &
     &         (b_cyl(1)*cos(this%alpha) - b_cyl(2)*sin(this%alpha))
 
         mse_get_modeled_signal(4) = mse_get_modeled_signal(2)                 &
     &                             / mse_get_modeled_signal(3)
 
         IF (btest(this%flags, mse_ratio_flag)) THEN
            mse_get_modeled_signal(1) = mse_get_modeled_signal(4)
         ELSE
            mse_get_modeled_signal(1) = atan2(mse_get_modeled_signal(2),       &
     &                                        mse_get_modeled_signal(3))
 
            IF (btest(this%flags, mse_degrees_flag)) THEN
               mse_get_modeled_signal(1) = mse_get_modeled_signal(1)           &
     &                                   / degree
            END IF
         END IF
 
         CALL this%scale_and_offset(a_model, mse_get_modeled_signal(1))
      ELSE
         mse_get_modeled_signal = last_value
      END IF
 
      CALL profiler_set_stop_time('mse_get_modeled_signal', start_time)
 
      END FUNCTION
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION mse_get_type(this)
      USE data_parameters
 
      IMPLICIT NONE
 
!  Declare Arguments
      CHARACTER (len=data_name_length) :: mse_get_type
      CLASS (mse_class), INTENT(in)    :: this
 
!  local variables
      REAL (rprec)                     :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      mse_get_type = 'mse'
 
      CALL profiler_set_stop_time('mse_get_type', start_time)
 
      END FUNCTION
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      SUBROUTINE mse_get_header(this, header)
      USE data_parameters
 
      IMPLICIT NONE
 
!  Declare Arguments
      class(mse_class), INTENT(in)                  :: this
      CHARACTER (len=data_name_length), DIMENSION(7), INTENT(inout) ::         &
     &   header
 
!  local variables
      REAL (rprec)                                   :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      header(1) = 'e_h'
      header(2) = 'e_v'
      header(3) = 'ratio'
 
      header(4) = 'model_sig(1)'
      header(5) = 'model_sig(2)'
      header(6) = 'model_sig(3)'
      header(7) = 'model_sig(4)'
 
      CALL profiler_set_stop_time('mse_get_header', start_time)
 
      END SUBROUTINE
 
!*******************************************************************************
!  PRIVATE
!*******************************************************************************
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION mse_angle_to_bt_proj(this, vec)
 
      IMPLICIT NONE
 
!  Declare Arguments
      REAL (rprec)                  :: mse_angle_to_bt_proj
      CLASS (mse_class), INTENT(in) :: this
      REAL (rprec), DIMENSION(3)    :: vec
 
!  local variables
      REAL (rprec), DIMENSION(2)    :: temp_vec
      REAL (rprec), DIMENSION(2)    :: norm_vec
      REAL (rprec)                  :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
!  Unit vector in the direction of
      temp_vec = (/ -this%point(2), this%point(1) /)                           &
     &         / sqrt(dot_product(this%point(1:2), this%point(1:2)))
 
!  Compute the angle.
      norm_vec = vec(1:2)/sqrt(dot_product(vec(1:2), vec(1:2)))
      mse_angle_to_bt_proj = acos(dot_product(temp_vec, norm_vec))
 
      CALL profiler_set_stop_time('mse_angle_to_bt_proj', start_time)
 
      END FUNCTION
 
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      FUNCTION mse_angle_to_horizontal(vec)
 
      IMPLICIT NONE
 
!  Declare Arguments
      REAL (rprec)               :: mse_angle_to_horizontal
      REAL (rprec), DIMENSION(3) :: vec
 
!  local variables
      REAL (rprec)               :: h
      REAL (rprec)               :: start_time
 
!  Start of executable code
      start_time = profiler_get_start_time()
 
      h = sqrt(dot_product(vec, vec))
      mse_angle_to_horizontal = asin(vec(3)/h)
      mse_angle_to_horizontal = pi/2.0 - mse_angle_to_horizontal
 
      CALL profiler_set_stop_time('mse_angle_to_horizontal', start_time)
 
      END FUNCTION
 
      END MODULE