text:
!>
!! @file
!! This file is part of SeisSol.
!!
!! @author Atanas Atanasov (atanasoa AT in.tum.de, http://www5.in.tum.de/wiki/index.php/Atanas_Atanasov)
!! @author Alice Gabriel (gabriel AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/gabriel)
!! @author Christian Pelties (pelties AT geophysik.uni-muenchen.de, http://www.geophysik.uni-muenchen.de/Members/pelties)
!! @author Sebastian Rettenberger (sebastian.rettenberger AT tum.de, http://www5.in.tum.de/wiki/index.php/Sebastian_Rettenberger)
!!
!! @section LICENSE
!! Copyright (c) 2013-2017, SeisSol Group
!! All rights reserved.
!!
!! Redistribution and use in source and binary forms, with or without
!! modification, are permitted provided that the following conditions are met:
!!
!! 1. Redistributions of source code must retain the above copyright notice,
!!    this list of conditions and the following disclaimer.
!!
!! 2. Redistributions in binary form must reproduce the above copyright notice,
!!    this list of conditions and the following disclaimer in the documentation
!!    and/or other materials provided with the distribution.
!!
!! 3. Neither the name of the copyright holder nor the names of its
!!    contributors may be used to endorse or promote products derived from this
!!    software without specific prior written permission.
!!
!! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
!! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
!! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
!! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
!! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!! POSSIBILITY OF SUCH DAMAGE.
!!
!! @section DESCRIPTION
!! Routines handling fault output

#include <Initializer/preProcessorMacros.fpp>

  MODULE faultoutput_mod
  !---------------------------------------------------------------------------!
  USE TypesDef
  !---------------------------------------------------------------------------!
  IMPLICIT NONE
  PRIVATE
  !---------------------------------------------------------------------------!
  ! Public procedures and functions
  !
  INTERFACE faultoutput
     MODULE PROCEDURE faultoutput
  END INTERFACE
  !
  !---------------------------------------------------------------------------!
  PUBLIC   :: faultoutput
  public   :: calc_FaultOutput
  PRIVATE  :: write_FaultOutput_elementwise
  PRIVATE  :: write_FaultOutput_atPickpoint
  !---------------------------------------------------------------------------!

CONTAINS

!> Fault output interface to manage output calculation and type of format
!<
  SUBROUTINE faultoutput(EQN, DISC, MESH, IO, MPI, MaterialVal, BND, time, dt)

      !-------------------------------------------------------------------------!
      USE JacobiNormal_mod
      USE magnitude_output_mod
      !-------------------------------------------------------------------------!
      IMPLICIT NONE
      !-------------------------------------------------------------------------!
      ! Argument list declaration
      TYPE(tEquations)         :: EQN
      TYPE(tDiscretization)    :: DISC
      TYPE(tUnstructMesh)      :: MESH
      TYPE(tInputOutput)       :: IO                                            ! IO structure
      TYPE(tMPI)               :: MPI                                           ! MPI
      TYPE(tBoundary)          :: BND                                           ! BND    data structure
      REAL                     :: MaterialVal(MESH%nElem,EQN%nBackgroundVar)    ! Local Mean Values
      REAL                     :: dt, time                                      ! Timestep and time
	  ! local variable declaration
      LOGICAL                  :: isOnPickpoint
      LOGICAL                  :: isOnElementwise
      LOGICAL                  :: isOnFaultoutput
      INTEGER                  :: i, rank_int                                   ! Loop variables
      !-------------------------------------------------------------------------!
      isOnPickpoint = .FALSE.
      isOnElementwise = .FALSE.

      !-------------------------------------------------------------------------!
      !
      ! Calculate Fault Output of the previous timestep
      ! Here, because the complete updated dgvar value of the (MPI-)Neighbor is needed
      ! Note that this causes a dt timeshift in the DR output routines
      !
      !
      IF (DISC%DynRup%energy_rate_output_on.EQ.1) THEN
         IF ( MOD(DISC%iterationstep,DISC%DynRup%energy_rate_printtimeinterval).EQ.0 &
         .OR. (DISC%EndTime-time).LE.(dt*1.005d0) ) THEN
            CALL energy_rate_output(MaterialVal,time,DISC,MESH,MPI,IO)
         ENDIF
      ENDIF

      SELECT CASE(DISC%DynRup%OutputPointType)
       ! For historical reasons fault output DISC%DynRup%OutputPointType= 3 or 4 or 5
       ! Case 0 means no fault output
       !CASE(0)
       !  logInfo(*) 'Fault output case (0).'
       !  RETURN

       ! Fault Receiver Output
       CASE(3)
         ! Check if output is desired
         IF (.NOT. DISC%DynRup%DynRup_out_atPickpoint%DR_pick_output ) THEN
             RETURN
         ENDIF
         !
         ! Check if output for this time step is desired
         ! (note that fault computation is always one behind due to MPI communication = -1!)
         IF ( MOD(DISC%iterationstep-1,DISC%DynRup%DynRup_out_atPickpoint%printtimeinterval).EQ.0 &
         .OR. (DISC%EndTime-time).LE.(dt*1.005d0) ) THEN
            IF (DISC%iterationstep.EQ.0) RETURN ! not the iteration 0
            CONTINUE
         ! print always first timestep
         ELSEIF(DISC%iterationstep .EQ. 1) THEN
            CONTINUE
         ELSE
            RETURN
         ENDIF
         CALL calc_FaultOutput(DISC%DynRup%DynRup_out_atPickpoint, DISC, EQN, MESH, MaterialVal, BND, time)
         CALL write_FaultOutput_atPickpoint(EQN, DISC, MESH, IO, MPI, MaterialVal, BND, time, dt)

       !Fault Element Output
       CASE(4)
         ! fault output iterations, including iteration 1 and last timestep
         IF (DISC%DynRup%DynRup_out_elementwise%printIntervalCriterion.EQ.1) THEN
            IF ( MOD(DISC%iterationstep-1,DISC%DynRup%DynRup_out_elementwise%printtimeinterval).EQ.0.OR. &
            ((min(DISC%EndTime,dt*DISC%MaxIteration)-time).LE.(dt*1.005d0))) THEN
               CONTINUE
!            ELSEIF(DISC%iterationstep .EQ. 0 ) THEN  ! uncomment to print setup at time 0 for intial stresses
!               CONTINUE
            ELSE  ! all timesteps without fault output
               RETURN
            ENDIF
         ELSE
            IF ((DISC%iterationstep.EQ.1).OR. ((DISC%iterationstep.GE.2) .AND. &
            ((MOD(time,DISC%DynRup%DynRup_out_elementwise%printtimeinterval_sec).LE.0.5*dt).OR. &
            ((ceiling(time/DISC%DynRup%DynRup_out_elementwise%printtimeinterval_sec)*DISC%DynRup%DynRup_out_elementwise%printtimeinterval_sec-time).LT.0.5*dt).OR. &
            ((min(DISC%EndTime,dt*DISC%MaxIteration)-time).LE.(dt*1.005d0))))) THEN
               CONTINUE
!            ELSEIF(DISC%iterationstep .EQ. 0 ) THEN  ! uncomment to print setup at time 0 for intial stresses
!               CONTINUE
            ELSE  ! all timesteps without fault output
               RETURN
            ENDIF
         ENDIF
         ! remember that fault output was written here
         isOnElementwise=.TRUE.
       ! combines option 3 and 4: output at individual stations and at complete fault
       CASE(5)
         ! check time for fault receiver output
         IF (.NOT. DISC%DynRup%DynRup_out_atPickpoint%DR_pick_output ) THEN
             CONTINUE
         ELSEIF ( MOD(DISC%iterationstep-1,DISC%DynRup%DynRup_out_atPickpoint%printtimeinterval).EQ.0 &
         .OR. (min(DISC%EndTime,dt*DISC%MaxIteration)-time).LE.(dt*1.005d0) ) THEN
            isOnPickpoint = .TRUE.
         ENDIF
         ! check time for fault plane output
         IF (DISC%DynRup%DynRup_out_elementwise%printIntervalCriterion.EQ.1) THEN
            IF ( MOD(DISC%iterationstep-1,DISC%DynRup%DynRup_out_elementwise%printtimeinterval).EQ.0.OR. &
            ((min(DISC%EndTime,dt*DISC%MaxIteration)-time).LE.(dt*1.005d0))) THEN
               isOnElementwise=.TRUE.
            ENDIF
         ELSE
            IF ((DISC%iterationstep.EQ.1).OR. ((DISC%iterationstep.GE.2) .AND. &
            ((MOD(time,DISC%DynRup%DynRup_out_elementwise%printtimeinterval_sec).LE.0.5*dt).OR. &
            ((ceiling(time/DISC%DynRup%DynRup_out_elementwise%printtimeinterval_sec)*DISC%DynRup%DynRup_out_elementwise%printtimeinterval_sec-time).LT.0.5*dt).OR. &
            ((min(DISC%EndTime,dt*DISC%MaxIteration)-time).LE.(dt*1.005d0))))) THEN
               isOnElementwise=.TRUE.
            ENDIF
         ENDIF

         ! print always first timestep
         IF(DISC%iterationstep .EQ. 1) THEN
           isOnPickpoint = .TRUE.
           isOnElementwise=.TRUE.
         ENDIF
         !
         IF (isOnPickpoint) THEN
           CALL calc_FaultOutput(DISC%DynRup%DynRup_out_atPickpoint, DISC, EQN, MESH, MaterialVal, BND, time)
           CALL write_FaultOutput_atPickpoint(EQN, DISC, MESH, IO, MPI, MaterialVal, BND, time, dt)
         ENDIF
         !
       CASE DEFAULT
          ! no output
          CONTINUE
      END SELECT
      !
  END SUBROUTINE faultoutput
!
!> Fault output calculation at specific positions (receiver and elementwise)
!<
  SUBROUTINE calc_FaultOutput( DynRup_output, DISC, EQN, MESH, MaterialVal, BND, time )
    use  f_ftoc_bind_interoperability
    USE create_fault_rotationmatrix_mod, only: create_strike_dip_unit_vectors
    use iso_c_binding, only: c_loc

    !-------------------------------------------------------------------------!
    USE common_operators_mod
    USE JacobiNormal_mod
    USE DGBasis_mod
    USE NucleationFunctions_mod
    !-------------------------------------------------------------------------!
    IMPLICIT NONE
    !-------------------------------------------------------------------------!
    ! Argument list declaration
    TYPE(tEquations)              :: EQN                                      !< EQN global variable
    TYPE(tDiscretization), TARGET :: DISC                                     !< DISC global variable
    TYPE(tUnstructMesh)           :: MESH                                     !< MESH global variable
    TYPE(tInputOutput)            :: IO                                       !< IO structure
    TYPE(tMPI)                    :: MPI                                      !< MPI
    TYPE(tBoundary)               :: BND                                      !< BND    data structure
    REAL                          :: MaterialVal(MESH%nElem,EQN%nBackgroundVar) !< Local Mean Values
    REAL                          :: time                                     !< time
    !-------------------------------------------------------------------------!
    ! Local variable declaration                                              !
    TYPE(tDynRup_output), target  :: DynRup_output                            !< Output data for Dynamic Rupture processes
    INTEGER :: iElem                                                          ! Element number                        !
    INTEGER :: iObject, MPIIndex, MPIIndex_DR
    INTEGER :: iDegFr                                                         ! Degree of freedom                     !
    INTEGER :: iVar                                                           ! Variable number                       !
    INTEGER :: iSide
    INTEGER :: iFace
    INTEGER :: iNeighbor, NeigBndGP
    INTEGER :: iLocalNeighborSide
    INTEGER :: NeighborPoly, LocPoly, MaxDegFr, iBndGP
    INTEGER :: NeighborDegFr, LocDegFr
    INTEGER :: OutVars, nOutPoints
    INTEGER :: i,j,m,k,iOutPoints                                             ! Loop variables                   !
    INTEGER :: SubElem, number_of_subtriangles                                ! elementwise fault refinement parameters
    REAL    :: NormalVect_n(3)                                                ! Normal vector components         !
    REAL    :: NormalVect_s(3)                                                ! Normal vector components         !
    REAL    :: NormalVect_t(3)                                                ! Normal vector components         !
    REAL    :: T(EQN%nVar,EQN%nVar)                                           ! Rotation matrix
    REAL    :: iT(EQN%nVar,EQN%nVar)                                          ! Rotation matrix
    REAL    :: V1(DISC%Galerkin%nDegFr,EQN%nVar)                              ! Reference state     (ref)
    REAL    :: V2(DISC%Galerkin%nDegFr,EQN%nVar)                              ! Reference state     (ref)
    REAL    :: Stress(6)                                                      ! The background stress tensor in vector form at a single BndGP
    REAL    :: SideVal(EQN%nVar), SideVal2(EQN%nVar)
    REAL    :: phi(2),rho, rho_neig, mu, mu_neig, lambda, lambda_neig
    REAL    :: LocXYStress, LocXZStress, TracXY, TracXZ, LocSRs, LocSRd
    REAL    :: w_speed(EQN%nNonZeroEV), TracEla, Trac, Strength, LocU, LocP, w_speed_neig(EQN%nNonZeroEV)
    REAL    :: S_0,P_0,S_XY,S_XZ,P_f, Temperature
    REAL    :: MuVal, cohesion, LocSV
    REAL    :: LocYY, LocZZ, LocYZ                                            ! temporary stress values
    REAL    :: tmp_mat(1:6), LocMat(1:6), TracMat(1:6)                        ! temporary stress tensors
    REAL    :: rotmat(1:6,1:6)                                                ! Rotation matrix
    REAL    :: TmpMat(EQN%nBackgroundVar)                                     ! temporary material values
    REAL    :: NorDivisor,ShearDivisor,UVelDivisor
    REAL    :: strike_vector(1:3), dip_vector(1:3), crossprod(1:3) !for rotation of Slip from local to strike, dip coordinate
    REAL    :: cos1, sin1, scalarprod
    REAL, PARAMETER    :: ZERO = 0.0D0
    ! Parameters used for calculating Vr
    LOGICAL :: compute_Vr
    INTEGER :: nDegFr2d, jBndGP, i1, j1
    REAL    :: chi, tau, Slowness, dt_dchi, dt_dtau, Vr
    REAL    :: dt_dx1, dt_dy1
    REAL    :: Tnuc, Gnuc, eta
    REAL    :: xV(4), yV(4), zV(4)
    REAL    :: xab(3), xac(3), grad2d(2,2), JacobiT2d(2,2)
    REAL    :: NucleationStressXYZ(1:6), NucleationStressLocalCS(1:6)
    REAL, ALLOCATABLE  :: projected_RT(:)
    real, allocatable  :: phiAtPoint(:)
    real, dimension( NUMBER_OF_BASIS_FUNCTIONS, NUMBER_OF_QUANTITIES ) :: DOFiElem_ptr ! no: it's not a pointer..
    real, dimension( NUMBER_OF_BASIS_FUNCTIONS, NUMBER_OF_QUANTITIES ) :: DOFiNeigh_ptr ! no pointer again
    !-------------------------------------------------------------------------!
    INTENT(IN)    :: BND, DISC, EQN, MESH, MaterialVal, time
    INTENT(INOUT) :: DynRup_output
    !-------------------------------------------------------------------------!
    !
    ! register epik/scorep function
    EPIK_FUNC_REG("calc_faultoutput")
    SCOREP_USER_FUNC_DEFINE()
    EPIK_FUNC_START()
    SCOREP_USER_FUNC_BEGIN("calc_faultoutput")
    !
    nOutPoints = DynRup_output%nDR_pick                                        ! number of output receivers for this MPI domain
    ! in case of refined fault output elements, the rotation matrices are only stored for every mother tet
    ! in pickpoint case SubElem = 1, in refined elementwise case SubElem = number of subtets per element
    IF (DISC%DynRup%OutputPointType.EQ.4) THEN
        IF(DISC%DynRup%DynRup_out_elementwise%refinement_strategy.eq.1) THEN
            number_of_subtriangles=3
        ELSE IF (DISC%DynRup%DynRup_out_elementwise%refinement_strategy.eq.2) THEN
            number_of_subtriangles=4
        ELSE
            number_of_subtriangles=1
        ENDIF
        SubElem = number_of_subtriangles**DISC%DynRup%DynRup_out_elementwise%refinement
    ELSE
        SubElem = 1
    ENDIF
    DO iOutPoints = 1,nOutPoints                                               ! loop over number of output receivers for this domain
          !
          iFace               = DynRup_output%RecPoint(iOutPoints)%index       ! current receiver location
          iElem               = MESH%Fault%Face(iFace,1,1)
          iSide               = MESH%Fault%Face(iFace,2,1)
          !
          iNeighbor           = MESH%Fault%Face(iFace,1,2)
          iLocalNeighborSide  = MESH%Fault%Face(iFace,2,2)
          !
          w_speed(:)      = DISC%Galerkin%WaveSpeed(iElem,:)
          rho             = MaterialVal(iElem,1)
          !
          if( iElem == 0 ) then
            call c_interoperability_getNeighborDofsFromDerivatives( i_meshId = iNeighbor, &
                                                                    i_faceId = iLocalNeighborSide, &
                                                                    o_dofs   = dofiElem_ptr )
          else
            call c_interoperability_getDofsFromDerivatives( i_meshId = iElem, &
                                                            o_dofs   = DOFiElem_ptr)
          endif

          IF (iNeighbor == 0) THEN
            ! iNeighbor is in the neighbor domain
            ! The neighbor element belongs to a different MPI domain
            iObject  = MESH%ELEM%BoundaryToObject(iSide,iElem)
            MPIIndex = MESH%ELEM%MPINumber(iSide,iElem)
            call c_interoperability_getNeighborDofsFromDerivatives( i_meshId = iElem, &
                                                                    i_faceId = iSide, &
                                                                    o_dofs   = DOFiNeigh_ptr )

            ! Bimaterial case only possible for elastic isotropic materials
            TmpMat(:)   = BND%ObjMPI(iObject)%NeighborBackground(:,MPIIndex)
            rho_neig    = TmpMat(1)
            mu_neig     = TmpMat(2)
            lambda_neig = TmpMat(3)
            w_speed_neig(1) = SQRT((lambda_neig+2.0D0*mu_neig)/rho_neig)  ! Will only work in elastic isotropic cases
            w_speed_neig(2) = SQRT(mu_neig/rho_neig)
            w_speed_neig(3) = w_speed_neig(2)
          ELSE
            ! normal case: iNeighbor present in local domain
            call c_interoperability_getDofsFromDerivatives( i_meshId = iNeighbor, &
                                                            o_dofs   = DOFiNeigh_ptr )
            w_speed_neig(:) = DISC%Galerkin%WaveSpeed(iNeighbor,:)
            rho_neig        = MaterialVal(iNeighbor,1)
          ENDIF
          !
          V1(:,:)=0.
          V2(:,:)=0.
          !
          ! currently no p-adaptivity
          ! IF(DISC%Galerkin%pAdaptivity.GE.1) THEN
          !   LocPoly  = INT(DISC%Galerkin%LocPoly(iElem))
          ! ELSE
          !   LocPoly  = DISC%Galerkin%nPoly
          ! ENDIF
          LocPoly  = DISC%Galerkin%nPoly
          LocDegFr = (LocPoly+1)*(LocPoly+2)*(LocPoly+3)/6.0D0
          !
          ! Local side's normal and tangential vectors
          NormalVect_n = MESH%Fault%geoNormals(1:3,iFace)
          NormalVect_s = MESH%Fault%geoTangent1(1:3,iFace)
          NormalVect_t = MESH%Fault%geoTangent2(1:3,iFace)
          !
          ! Rotate DoF
          CALL RotationMatrix3D(NormalVect_n,NormalVect_s,NormalVect_t,T(:,:),iT(:,:),EQN)
          DO iDegFr=1,LocDegFr
            V1(iDegFr,:)=MATMUL(iT(:,:),dofiElem_ptr(iDegFr,1:EQN%nVar))
            V2(iDegFr,:)=MATMUL(iT(:,:),DOFiNeigh_ptr(iDegFr,1:EQN%nVar))
          ENDDO
          !
          ! load nearest boundary GP: iBndGP
          iBndGP = DynRup_output%OutInt(iOutPoints,1)

          !Background stress rotation to face's reference system
          Stress = EQN%InitialStressInFaultCS(iBndGP,:,iFace)
          !
          MuVal = DISC%DynRup%output_Mu(iBndGP,iFace)
          LocSV = DISC%DynRup%output_StateVar(iBndGP,iFace) ! load state variable of RS for output
          cohesion  = DISC%DynRup%cohesion(iBndGP,iFace)
          S_XY  = Stress(4)
          S_XZ  = Stress(6)
          P_0   = Stress(1)
          IF (DISC%DynRup%ThermalPress.EQ.1) THEN
             P_f = DISC%DynRup%TP(iBndGP, iFace, 2)
             Temperature = DISC%DynRup%TP(iBndGP, iFace, 1)
          ELSE
             P_f = 0.0
          ENDIF

          if (EQN%FL.eq.33) then
             !case of ImposedSlipRateOnDRBoundary 'friction law': we add the additional stress to the fault output
             !to show the imposed SR
             eta = (w_speed(2)*rho*w_speed_neig(2)*rho_neig) / (w_speed(2)*rho + w_speed_neig(2)*rho_neig)
             Gnuc = Calc_SmoothStep(time, DISC%DynRup%t_0)
             S_XY = S_XY - eta * EQN%NucleationStressInFaultCS(iBndGP,1,iFace)*Gnuc
             S_XZ = S_XZ - eta * EQN%NucleationStressInFaultCS(iBndGP,2,iFace)*Gnuc
          endif

          !
          ! Obtain values at output points
          SideVal  = 0.
          SideVal2 = 0.
          DO iDegFr = 1, LocDegFr
             ! Basis functions for the automatic evaluation of DOFs at fault output nodes
             phi(:) = DynRup_output%OutEval(iOutPoints,1,iDegFr,:)
             SideVal(:)  = SideVal(:)  + V1(iDegFr,:)*phi(1)
             SideVal2(:) = SideVal2(:) + V2(iDegFr,:)*phi(2)
          ENDDO
          !
          ! Compute divisors
          NorDivisor   = 1.0D0/(w_speed_neig(1)*rho_neig+w_speed(1)*rho)
          ShearDivisor = 1.0D0/(w_speed_neig(2)*rho_neig+w_speed(2)*rho)
          UVelDivisor  = 1.0D0/(w_speed(1)*rho)
          !
          ! Traction (elastic response, not at sliding fault!)
          LocXYStress     =  SideVal(4) + &
                            (((SideVal2(4)-SideVal(4)) + &
                            w_speed_neig(2)*rho_neig*(SideVal2(8)-SideVal(8)))* &
                            w_speed(2)*rho ) * ShearDivisor
          LocXZStress     =  SideVal(6) + &
                            (((SideVal2(6)-SideVal(6)) + &
                            w_speed_neig(2)*rho_neig*(SideVal2(9)-SideVal(9)))* &
                            w_speed(2)*rho ) * ShearDivisor
          ! Normal stress
          LocP = SideVal(1)+(((SideVal2(1)-SideVal(1))+w_speed_neig(1)*rho_neig*(SideVal2(7)-SideVal(7)))* &
                 w_speed(1)*rho) * NorDivisor
          ! Normal velocity
          LocU = SideVal(7)+ (LocP-SideVal(1)) * UVelDivisor

          ! compute missing sigma values for rotation
          LocYY = SideVal(2) + (LocP-SideVal(1)) * (1.0D0-2.0D0*w_speed(2)**2/w_speed(1)**2)
          LocZZ = SideVal(3) + (LocP-SideVal(1)) * (1.0D0-2.0D0*w_speed(2)**2/w_speed(1)**2)
          LocYZ = SideVal(5)

          ! absolute shear stress
          TracEla = SQRT((S_XY + LocXYStress)**2 + (S_XZ + LocXZStress)**2)

          ! compute fault strength
          SELECT CASE(EQN%FL)
          CASE DEFAULT
            ! linear slip weakening
            Strength = -MuVal*MIN(LocP+P_0 -P_f,ZERO) - cohesion
          CASE(3,4)
             ! rate and state (once everything is tested and cohesion works for RS, this option could be merged to default)
             Strength = -MuVal*MIN(LocP+P_0 -P_f, ZERO)
          CASE(6)
            ! exception for bimaterial with LSW case
            ! modify strength according to prakash clifton
            ! in the first step we just take the Strength values of the nearest GP
            Strength = DISC%DynRup%output_Strength(iBndGP,iFace)
          END SELECT
          !
          ! Coulomb's model for obtaining actual traction
          IF(TracEla.GT.ABS(Strength)) THEN
              TracXY = ((S_XY+LocXYStress)/TracEla)*Strength
              TracXZ = ((S_XZ+LocXZStress)/TracEla)*Strength
              !
              ! update stress change
              TracXY = TracXY - S_XY
              TracXZ = TracXZ - S_XZ
          ELSE
              Trac = TracEla
              TracXY = LocXYStress
              TracXZ = LocXZStress
          ENDIF
          !
          ! Rotate back to xyz-system
          ! create dummy matrix
          ! tmp_mat in xyz system:
          tmp_mat    = 0.0D0
          tmp_mat(1) = LocP
          tmp_mat(2) = LocYY
          tmp_mat(3) = LocZZ
          tmp_mat(4) = TracXY ! modified traction !
          tmp_mat(5) = LocYZ
          tmp_mat(6) = TracXZ ! modified traction !
          tmp_mat    = MATMUL(T(1:6,1:6),tmp_mat)

          ! rotate into fault system
          rotmat  = DynRup_output%rotmat((iOutPoints-1)/(SubElem)+1,1:6,1:6) ! in pickpoint case SubElem = 1, in refined elementwise case SubElem = number of subtets per element
          TracMat = MATMUL(rotmat,tmp_mat)

          ! Do it again but assume no slip occured
          tmp_mat    = 0.0D0
          tmp_mat(1) = LocP
          tmp_mat(2) = LocYY
          tmp_mat(3) = LocZZ
          tmp_mat(4) = LocXYStress ! traction of standard wave propagation without slip!
          tmp_mat(5) = LocYZ
          tmp_mat(6) = LocXZStress ! traction of standard wave propagation without slip!
          tmp_mat    = MATMUL(T(1:6,1:6),tmp_mat) ! tmp_mat in xyz system

          ! rotate into fault system
          LocMat = MATMUL(rotmat,tmp_mat)

          CALL create_strike_dip_unit_vectors(NormalVect_n, strike_vector, dip_vector)

          ! sliprate
          if (DISC%DynRup%SlipRateOutputType .eq. 1) then
            LocSRs = -(1.0D0/(w_speed(2)*rho)+1.0D0/(w_speed_neig(2)*rho_neig))*(TracMat(4)-LocMat(4))
            LocSRd = -(1.0D0/(w_speed(2)*rho)+1.0D0/(w_speed_neig(2)*rho_neig))*(TracMat(6)-LocMat(6))
          else
            LocSRs = dot_product(SideVal2(8) * NormalVect_s + SideVal2(9) * NormalVect_t, strike_vector) - dot_product(SideVal(8) * NormalVect_s + SideVal(9) * NormalVect_t, strike_vector)
            LocSRd = dot_product(SideVal2(8) * NormalVect_s + SideVal2(9) * NormalVect_t, dip_vector   ) - dot_product(SideVal(8) * NormalVect_s + SideVal(9) * NormalVect_t, dip_vector   )
          end if

          ! TU 06.2015: average stress over the element in the considered direction
          ! (uncomment for rough geometry)
          !SideVal(1)  =  V1(1,1)*DynRup_output%OutEval(iOutPoints,1,1,1)
          !SideVal2(1) =  V2(1,1)*DynRup_output%OutEval(iOutPoints,1,1,2)
          !LocP = SideVal(1)+(((SideVal2(1)-SideVal(1))+w_speed_neig(1)*rho_neig*(SideVal2(7)-SideVal(7)))* &
          !       w_speed(1)*rho) * NorDivisor
          ! Store Values into Output vector OutVal

          if (EQN%FL.eq.33) then
             !case of ImposedSlipRateOnDRBoundary 'friction law': we plot the Stress from Godunov state, because we want to see the traction change from the imposed slip distribution
             TracMat(4)=LocMat(4)
             TracMat(6)=LocMat(6)
          endif

          OutVars = 0
          IF (DynRup_output%OutputMask(1).EQ.1) THEN
              OutVars = OutVars + 1
              DynRup_output%OutVal(iOutPoints,1,OutVars) = LocSRs !OutVars =1
              OutVars = OutVars + 1
              DynRup_output%OutVal(iOutPoints,1,OutVars) = LocSRd !OutVars =2
          ENDIF
          IF (DynRup_output%OutputMask(2).EQ.1) THEN
              OutVars = OutVars + 1
              DynRup_output%OutVal(iOutPoints,1,OutVars) = TracMat(4) !OutVars =3
              OutVars = OutVars + 1
              DynRup_output%OutVal(iOutPoints,1,OutVars) = TracMat(6) !OutVars =4
              OutVars = OutVars + 1
              DynRup_output%OutVal(iOutPoints,1,OutVars) = LocP-P_f !OutVars =5
          ENDIF
          IF (DynRup_output%OutputMask(3).EQ.1) THEN
              OutVars = OutVars + 1
              DynRup_output%OutVal(iOutPoints,1,OutVars) = LocU !OutVars =6
          ENDIF
          IF (DynRup_output%OutputMask(4).EQ.1) THEN
              OutVars = OutVars + 1
              DynRup_output%OutVal(iOutPoints,1,OutVars) = MuVal !OutVars =7
              OutVars = OutVars + 1
              DynRup_output%OutVal(iOutPoints,1,OutVars) = LocSV !OutVars =8
          ENDIF
          IF (DynRup_output%OutputMask(5).EQ.1) THEN

              !add transient nucleation stress to fault output
              IF ((EQN%FL.EQ.103) .OR. (EQN%FL.EQ.3) .OR. (EQN%FL.EQ.4)) THEN
                 NucleationStressXYZ = MATMUL(T(1:6,1:6), EQN%NucleationStressInFaultCS(iBndGP,:,iFace))
                 Tnuc = DISC%DynRup%t_0
                 Gnuc = Calc_SmoothStep(time, Tnuc)
                 NucleationStressLocalCS = MATMUL(rotmat, NucleationStressXYZ)*Gnuc
              ELSE
                 NucleationStressLocalCS(:) = 0d0
              ENDIF

              OutVars = OutVars + 1
              DynRup_output%OutVal(iOutPoints,1,OutVars) = TracMat(4)+DISC%DynRup%DynRup_Constants(iOutPoints)%ts0 + NucleationStressLocalCS(4) !OutVars =9
              OutVars = OutVars + 1
              DynRup_output%OutVal(iOutPoints,1,OutVars) = TracMat(6)+DISC%DynRup%DynRup_Constants(iOutPoints)%td0 + NucleationStressLocalCS(6) !OutVars =10
              OutVars = OutVars + 1
              DynRup_output%OutVal(iOutPoints,1,OutVars) = LocP+DISC%DynRup%DynRup_Constants(iOutPoints)%p0 - P_f + NucleationStressLocalCS(1)!OutVars =11
          ENDIF

              IF (DynRup_output%OutputMask(6).EQ.1) THEN
                  ! TU 07.15 rotate Slip from face reference coordinate to (strike,dip, normal) reference cordinate
                  cos1 = dot_product(strike_vector(:),NormalVect_s(:))
                  crossprod(:) = strike_vector(:) .x. NormalVect_s(:)

                  scalarprod = dot_product(crossprod(:),NormalVect_n(:))
                  !TU 2.11.15 :cos1**2 can be greater than 1 because of rounding errors -> min
                  IF (scalarprod.GT.0) THEN
                      sin1=sqrt(1-min(1.0,cos1**2))
                  ELSE
                      sin1=-sqrt(1-min(1.0,cos1**2))
                  ENDIF

                  OutVars = OutVars + 1
                  DynRup_output%OutVal(iOutPoints,1,OutVars)  = cos1 * DISC%DynRup%output_Slip1(iBndGP,iFace) - sin1* DISC%DynRup%output_Slip2(iBndGP,iFace)
                  OutVars = OutVars + 1
                  DynRup_output%OutVal(iOutPoints,1,OutVars)  = sin1 * DISC%DynRup%output_Slip1(iBndGP,iFace) + cos1 * DISC%DynRup%output_Slip2(iBndGP,iFace)
              ENDIF

              IF (DynRup_output%OutputMask(7).EQ.1) THEN
                ! Thomas Ulrich 3.08.2015
                ! Calculation of rupture velocity from rupture time array
                ! using the 2d basis functions and their derivatives

                ! To compute accurately dt/dxi, all rupture_time have to be defined on the element
                compute_VR = .TRUE.
                DO jBndGP = 1,DISC%Galerkin%nBndGP
                   IF (DISC%DynRup%output_rupture_time(jBndGP,iFace).EQ.0) THEN
                      compute_VR = .FALSE.
                      Vr = 0d0
                   ENDIF
                ENDDO

                IF (compute_VR.EQ..TRUE.) THEN

                ! The rupture velocity is based on the gradient of the rupture time
                ! The gradient is not expected to be accurate for the side gauss points
                !! We consider then the next internal gauss point if the closest gauss point is external
                i1 = int((iBndGP-1)/(DISC%Galerkin%nPoly + 2))+1
                j1 = modulo((iBndGP-1),(DISC%Galerkin%nPoly + 2))+1
                IF (i1.EQ.1) THEN
                   i1=i1+1
                ELSE IF (i1.EQ.(DISC%Galerkin%nPoly + 2)) THEN
                   i1=i1-1
                ENDIF
                IF (j1.EQ.1) THEN
                   j1=j1+1
                ELSE IF (j1.EQ.(DISC%Galerkin%nPoly + 2)) THEN
                   j1=j1-1
                ENDIF
                iBndGP = (i1-1)*(DISC%Galerkin%nPoly + 2)+j1

                !Uncomment for a selecting only the central GP
                !i1  = int((DISC%Galerkin%nPoly + 2 -1)/2)+1
                !iBndGP = (i1-1)*(DISC%Galerkin%nPoly + 2)+i1

                !project rupture_time onto the basis functions
                nDegFr2d = (DISC%Galerkin%nPoly + 1)*(DISC%Galerkin%nPoly + 2)/2
                !projection of the rupture time on the basis functions
                ALLOCATE(projected_RT(nDegFr2d))
                ALLOCATE(phiAtPoint(2 * nDegFr2d))
                projected_RT(:)=0.

                DO jBndGP = 1,DISC%Galerkin%nBndGP
                  chi  = MESH%ELEM%BndGP_Tri(1,jBndGP)
                  tau  = MESH%ELEM%BndGP_Tri(2,jBndGP)
                  call c_interoperability_TriDubinerP(phiAtPoint, chi, tau, DISC%Galerkin%nPoly)
                  DO iDegFr = 1, nDegFr2d
                     projected_RT(iDegFr) = projected_RT(iDegFr) + &
                        DISC%Galerkin%BndGaussW_Tet(jBndGP) * &
                        DISC%DynRup%output_rupture_time(jBndGP,iFace) * phiAtPoint(iDegFr)
                  ENDDO
                ENDDO
                DO iDegFr = 1,nDegFr2d
                   projected_RT(iDegFr) = projected_RT(iDegFr) * &
                       c_interoperability_M2invDiagonal(iDegFr-1)
                ENDDO

                !calculation of the spatial derivatives of the rupture time
                chi  = MESH%ELEM%BndGP_Tri(1,iBndGP)
                tau  = MESH%ELEM%BndGP_Tri(2,iBndGP)
                call c_interoperability_gradTriDubinerP(phiAtPoint, chi, tau, DISC%Galerkin%nPoly)
                dt_dchi=0d0
                dt_dtau=0d0
                DO iDegFr = 1, nDegFr2d
                   dt_dchi = dt_dchi + projected_RT(iDegFr) * phiAtPoint(1 + 2*(iDegFr-1))
                   dt_dtau = dt_dtau + projected_RT(iDegFr) * phiAtPoint(2 + 2*(iDegFr-1))
                ENDDO

                !For calculating the 2d Jacobian matrix we will need the nodal coordinates
                IF (MESH%Fault%Face(iFace,1,1) == 0) THEN
                ! iElem is in the neighbor domain
                ! The neighbor element belongs to a different MPI domain
                iNeighbor           = MESH%Fault%Face(iFace,1,2)          ! iNeighbor denotes "-" side
                iLocalNeighborSide  = MESH%Fault%Face(iFace,2,2)
                iObject  = MESH%ELEM%BoundaryToObject(iLocalNeighborSide,iNeighbor)
                MPIIndex = MESH%ELEM%MPINumber(iLocalNeighborSide,iNeighbor)
                !
                xV(1:4) = BND%ObjMPI(iObject)%NeighborCoords(1,1:4,MPIIndex)
                yV(1:4) = BND%ObjMPI(iObject)%NeighborCoords(2,1:4,MPIIndex)
                zV(1:4) = BND%ObjMPI(iObject)%NeighborCoords(3,1:4,MPIIndex)
                ELSE
                !
                ! get vertices
                xV(1:4) = MESH%VRTX%xyNode(1,MESH%ELEM%Vertex(1:4,iElem))
                yV(1:4) = MESH%VRTX%xyNode(2,MESH%ELEM%Vertex(1:4,iElem))
                zV(1:4) = MESH%VRTX%xyNode(3,MESH%ELEM%Vertex(1:4,iElem))
                ENDIF

                SELECT CASE(iSide)
                 CASE(1)
                    xab = (/ xV(3)-xV(1), yV(3)-yV(1), zV(3)-zV(1) /)
                    xac = (/ xV(2)-xV(1), yV(2)-yV(1), zV(2)-zV(1) /)
                 CASE(2)
                    xab = (/ xV(2)-xV(1), yV(2)-yV(1), zV(2)-zV(1) /)
                    xac = (/ xV(4)-xV(1), yV(4)-yV(1), zV(4)-zV(1) /)
                 CASE(3)
                    xab = (/ xV(4)-xV(1), yV(4)-yV(1), zV(4)-zV(1) /)
                    xac = (/ xV(3)-xV(1), yV(3)-yV(1), zV(3)-zV(1) /)
                 CASE(4)
                    xab = (/ xV(3)-xV(2), yV(3)-yV(2), zV(3)-zV(2) /)
                    xac = (/ xV(4)-xV(2), yV(4)-yV(2), zV(4)-zV(2) /)
                ENDSELECT

                grad2d(1,1) = dot_product(NormalVect_s(:),xab)
                grad2d(1,2) = dot_product(NormalVect_t(:),xab)
                grad2d(2,1) = dot_product(NormalVect_s(:),xac)
                grad2d(2,2) = dot_product(NormalVect_t(:),xac)
                CALL MatrixInverse2x2(JacobiT2d,grad2d)

                dt_dx1 = JacobiT2d(1,1) * dt_dchi + JacobiT2d(1,2) * dt_dtau
                dt_dy1 = JacobiT2d(2,1) * dt_dchi + JacobiT2d(2,2) * dt_dtau
                Slowness = sqrt(dt_dx1**2 + dt_dy1**2)

                IF (Slowness.EQ.0d0) THEN
                  Vr=0d0
                ElSE
                  Vr=1d0/Slowness
                ENDIF
                DEALLOCATE(projected_RT)
                DEALLOCATE(phiAtPoint)
                iBndGP = DynRup_output%OutInt(iOutPoints,1)
                ENDIF

                OutVars = OutVars + 1
                DynRup_output%OutVal(iOutPoints,1,OutVars)  = Vr

              ENDIF
              IF (DynRup_output%OutputMask(8).EQ.1) THEN
                  OutVars = OutVars + 1
                  DynRup_output%OutVal(iOutPoints,1,OutVars) = DISC%DynRup%output_Slip(iBndGP,iFace)
              ENDIF
              IF (DynRup_output%OutputMask(9).EQ.1) THEN
                  OutVars = OutVars + 1
                  DynRup_output%OutVal(iOutPoints,1,OutVars) = DISC%DynRup%output_PeakSR(iBndGP,iFace)
              ENDIF
              IF (DynRup_output%OutputMask(10).EQ.1) THEN
                  OutVars = OutVars + 1
                  DynRup_output%OutVal(iOutPoints,1,OutVars) = DISC%DynRup%output_rupture_time(iBndGP,iFace)
              ENDIF
              IF (DynRup_output%OutputMask(11).EQ.1) THEN
                  OutVars = OutVars + 1
                  DynRup_output%OutVal(iOutPoints,1,OutVars) = DISC%DynRup%output_dynStress_time(iBndGP,iFace)
              ENDIF
              IF (DISC%DynRup%ThermalPress.EQ.1) THEN
                  IF (DynRup_output%OutputMask(12).EQ.1) THEN
                      OutVars = OutVars + 1
                      DynRup_output%OutVal(iOutPoints,1,OutVars) = DISC%DynRup%TP(iBndGP,iFace,2) !OutVars = 19
                      OutVars = OutVars + 1
                      DynRup_output%OutVal(iOutPoints,1,OutVars) = DISC%DynRup%TP(iBndGP,iFace,1) !OutVars = 20
                  ENDIF
              ENDIF

          ! Store output
          IF (DISC%DynRup%OutputPointType.NE.4) THEN
          DynRup_output%CurrentPick(iOutPoints) = DynRup_output%CurrentPick(iOutPoints) +1
          !
          DynRup_output%TmpTime(DynRup_output%CurrentPick(iOutPoints)) = time
          !
          ELSE
          DynRup_output%CurrentPick(iOutPoints) = 1
          ENDIF
          DynRup_output%TmpState(iOutPoints,DynRup_output%CurrentPick(iOutPoints),:) = DynRup_output%OutVal(iOutPoints,1,:)
          !
    ENDDO ! iOutPoints = 1,nOutPoints
    !
    CONTINUE

    EPIK_FUNC_END()
    SCOREP_USER_FUNC_END()
    !
  END SUBROUTINE calc_FaultOutput

  !
  !> Subroutine initializing the fault output
  !<
  SUBROUTINE write_FaultOutput_elementwise(EQN, DISC, MESH, IO, MPI, MaterialVal, BND, time, dt)
      !-------------------------------------------------------------------------!
      USE JacobiNormal_mod
      use FaultWriter
      use, intrinsic :: iso_c_binding
      !-------------------------------------------------------------------------!
      IMPLICIT NONE
      !-------------------------------------------------------------------------!
      ! Argument list declaration
      TYPE(tEquations)         :: EQN
      TYPE(tDiscretization)    :: DISC
      TYPE(tUnstructMesh)      :: MESH
      TYPE(tInputOutput)       :: IO                                            ! IO structure
      TYPE(tMPI)               :: MPI                                           ! MPI
      TYPE(tBoundary)          :: BND                                           ! BND    data structure
      REAL                     :: MaterialVal(MESH%nElem,EQN%nBackgroundVar)    ! Local Mean Values
      REAL                     :: dt, time
      !-------------------------------------------------------------------------!
      call writeFault(time)
  END SUBROUTINE ! write_FaultOutput_elementwise
  !
  !> writing fault output at pickpoints to files
  !<
  SUBROUTINE write_FaultOutput_atPickpoint(EQN, DISC, MESH, IO, MPI, MaterialVal, BND, time, dt)
    !-------------------------------------------------------------------------!
    USE JacobiNormal_mod
    !-------------------------------------------------------------------------!
    IMPLICIT NONE
    !-------------------------------------------------------------------------!
    ! Argument list declaration
    TYPE(tEquations)              :: EQN
    TYPE(tDiscretization)         :: DISC
    TYPE(tUnstructMesh)           :: MESH
    TYPE(tInputOutput)            :: IO                                         ! IO structure
    TYPE(tMPI)                    :: MPI                                        ! MPI
    TYPE(tBoundary)               :: BND                                        ! BND    data structure
    REAL                          :: MaterialVal(MESH%nElem,EQN%nBackgroundVar) ! Local Mean Values
    !-------------------------------------------------------------------------!
    ! Local variable declaration                                              !
    REAL    :: dt, time                                                       ! Timestep and time
    INTEGER :: stat
    INTEGER :: nOutPoints
    INTEGER :: iOutPoints, k                                                  ! Loop variables                   !
    REAL, ALLOCATABLE       :: TmpStateIO(:)                                  ! local variable to avoid fortran runtime copy
    CHARACTER (LEN=5)       :: cmyrank
    CHARACTER (len=200)     :: ptsoutfile
    !-------------------------------------------------------------------------!
    INTENT(IN)    :: EQN, MESH, MaterialVal, time
    INTENT(INOUT) :: DISC
    !-------------------------------------------------------------------------!
    !

    ALLOCATE ( TmpStateIO( UBOUND(DISC%DynRup%DynRup_out_atPickpoint%TmpState, 3) ) )
    nOutPoints = DISC%DynRup%DynRup_out_atPickpoint%nDR_pick                   ! number of output receivers for this MPI domain
    DO iOutPoints = 1,nOutPoints                                               ! loop over number of output receivers for this domain
          !
          IF(DISC%DynRup%DynRup_out_atPickpoint%CurrentPick(iOutPoints).GE.DISC%DynRup%DynRup_out_atPickpoint%MaxPickStore.OR.ABS(DISC%EndTime-time).LE.(dt*1.005d0)) THEN
#ifdef PARALLEL
            WRITE(cmyrank,'(I5.5)') MPI%myrank                                   ! myrank -> cmyrank
            WRITE(ptsoutfile, '(a,a15,i5.5,a1,a5,a4)') TRIM(IO%OutputFile),'-faultreceiver-',DISC%DynRup%DynRup_out_atPickpoint%RecPoint(iOutPoints)%globalreceiverindex,'-',TRIM(cmyrank),'.dat'!
#else
            WRITE(ptsoutfile, '(a,a15,i5.5,a4)') TRIM(IO%OutputFile),'-faultreceiver-',DISC%DynRup%DynRup_out_atPickpoint%RecPoint(iOutPoints)%globalreceiverindex,'.dat'
#endif
            !
            OPEN(UNIT     = DISC%DynRup%DynRup_out_atPickpoint%VFile(iOutPoints)                    , & !
                 FILE     = ptsoutfile                                       , & !
                 FORM     = 'FORMATTED'                                      , & !
                 STATUS   = 'OLD'                                            , & !
                 POSITION = 'APPEND'                                         , & !
                 RECL     = 80000                                            , & !
                 IOSTAT = stat                                                 ) !
            IF( stat.NE.0) THEN
               logError(*) 'cannot open ',ptsoutfile
               logError(*) 'Error status: ', stat
               call MPI_ABORT(MPI%commWorld, 134)
            END IF
            !
            DO k=1,DISC%DynRup%DynRup_out_atPickpoint%CurrentPick(iOutPoints)
               TmpStateIO(:) = DISC%DynRup%DynRup_out_atPickpoint%TmpState(iOutPoints,k,:)
               WRITE(DISC%DynRup%DynRup_out_atPickpoint%VFile(iOutPoints),*) DISC%DynRup%DynRup_out_atPickpoint%TmpTime(k), TmpStateIO(:)
            ENDDO
            !
            DISC%DynRup%DynRup_out_atPickpoint%CurrentPick(iOutPoints) = 0
            !
            CLOSE( DISC%DynRup%DynRup_out_atPickpoint%VFile(iOutPoints) ) ! not the same numbers as in elementwise output
            !
          ENDIF
    ENDDO ! iOutPoints = 1,DISC%DynRup%nOutPoints
    DEALLOCATE ( TmpStateIO )
  END SUBROUTINE write_FaultOutput_atPickpoint
 !
 END MODULE faultoutput_mod

max_stars_repo_name: LukasvdWiel/SeisSol

max_issues_repo_name: LukasvdWiel/SeisSol

max_forks_repo_name: LukasvdWiel/SeisSol


!******************************************************************************
!******************************************************************************
MODULE estimate_fortran

    !/* external modules    */

    USE optimizers_interfaces

    USE recording_estimation

    USE recording_warning

    USE shared_interface

    USE evaluate_fortran

    USE solve_fortran

#if MPI_AVAILABLE

    USE parallelism_constants

    USE parallelism_auxiliary

    USE mpi

#endif

    !/* setup   */

    IMPLICIT NONE

    PUBLIC

CONTAINS
!******************************************************************************
!******************************************************************************
SUBROUTINE fort_estimate(crit_val, success, message, optim_paras, optimizer_used, maxfun, num_procs, edu_spec, precond_spec, optimizer_options, num_types)

    !/* external objects    */

    TYPE(OPTIMIZER_COLLECTION), INTENT(INOUT) :: optimizer_options

    REAL(our_dble), INTENT(OUT)         :: crit_val

    TYPE(OPTIMPARAS_DICT), INTENT(IN)   :: optim_paras
    TYPE(EDU_DICT), INTENT(IN)          :: edu_spec

    TYPE(PRECOND_DICT), INTENT(IN)      :: precond_spec

    INTEGER(our_int), INTENT(IN)        :: num_procs
    INTEGER(our_int), INTENT(IN)        :: num_types
    INTEGER(our_int), INTENT(IN)        :: maxfun

    CHARACTER(225), INTENT(IN)          :: optimizer_used
    CHARACTER(150), INTENT(OUT)         :: message

    LOGICAL, INTENT(OUT)                :: success

    !/* internal objects    */

    REAL(our_dble)                  :: x_optim_free_scaled_start(num_free)

    REAL(our_dble)                  :: x_optim_free_unscaled_start(num_free)

    INTEGER(our_int)                :: iter

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

    IF (num_procs == 1) THEN
        criterion_function => fort_criterion_scalar
    ELSE
        criterion_function => fort_criterion_parallel
    END IF

    ! We need to determine the number of observations we have for each individual before any evaluation of the criterion function can be conducted. However, we only need to do this here if we are not letting the slaves do all the work. In that case the dataset is not read on the head node.
    IF (num_procs == 1) num_obs_agent = get_num_obs_agent(data_est)

    ! Some ingredients for the evaluation of the criterion function need to be created once and shared globally.
    CALL get_optim_paras(x_all_start, optim_paras, .True.)

    CALL fort_create_state_space(states_all, states_number_period, mapping_state_idx, num_periods, num_types, edu_spec)

    CALL get_optim_paras(x_optim_free_unscaled_start, optim_paras, .False.)

    CALL get_precondition_matrix(precond_matrix, precond_spec, maxfun, x_optim_free_unscaled_start)

    x_optim_free_scaled_start = apply_scaling(x_optim_free_unscaled_start, precond_matrix, 'do')
    x_optim_bounds_free_scaled(1, :) = apply_scaling(x_optim_bounds_free_unscaled(1, :), precond_matrix, 'do')
    x_optim_bounds_free_scaled(2, :) = apply_scaling(x_optim_bounds_free_unscaled(2, :), precond_matrix, 'do')

    CALL record_estimation(precond_matrix, x_optim_free_unscaled_start, optim_paras, .False.)

    CALL auto_adjustment_optimizers(optimizer_options, optimizer_used)

    crit_estimation = .True.

    IF (maxfun == zero_int) THEN

        success = .True.
        message = 'Single evaluation of criterion function at starting values.'

        CALL record_estimation('Start')
        crit_val = criterion_function(x_optim_free_scaled_start)
        CALL record_estimation('Finish')

    ELSEIF (optimizer_used == 'FORT-NEWUOA') THEN

        CALL newuoa(criterion_function, x_optim_free_scaled_start, optimizer_options%newuoa%npt, optimizer_options%newuoa%rhobeg, optimizer_options%newuoa%rhoend, zero_int, MIN(maxfun, optimizer_options%newuoa%maxfun), success, message)

    ELSEIF (optimizer_used == 'FORT-BOBYQA') THEN

        ! The BOBYQA algorithm might adjust the starting values. So we simply make sure that the very first evaluation of the criterion function is at the actual starting values.
        crit_val = criterion_function(x_optim_free_scaled_start)
        CALL bobyqa(criterion_function, x_optim_free_scaled_start, optimizer_options%bobyqa%npt, optimizer_options%bobyqa%rhobeg, optimizer_options%bobyqa%rhoend, zero_int, MIN(maxfun, optimizer_options%bobyqa%maxfun), success, message)

    ELSEIF (optimizer_used == 'FORT-BFGS') THEN
        dfunc_eps = optimizer_options%bfgs%eps
        CALL dfpmin(criterion_function, fort_dcriterion, x_optim_free_scaled_start, optimizer_options%bfgs%gtol, optimizer_options%bfgs%maxiter, optimizer_options%bfgs%stpmx, maxfun, success, message, iter)
        dfunc_eps = -HUGE_FLOAT
    END IF

    crit_estimation = .False.

    CALL record_estimation(success, message)

    CALL record_estimation()

END SUBROUTINE
!******************************************************************************
!******************************************************************************
FUNCTION fort_criterion_scalar(x_optim_free_scaled)

    !/* external objects    */

    REAL(our_dble), INTENT(IN)      :: x_optim_free_scaled(:)
    REAL(our_dble)                  :: fort_criterion_scalar

    !/* internal objects    */

    REAL(our_dble)                  :: x_optim_all_unscaled(num_paras)
    REAL(our_dble)                  :: x_optim_free_unscaled(num_free)
    REAL(our_dble)                  :: contribs(num_agents_est)
    REAL(our_dble)                  :: start

    INTEGER(our_int)                :: dist_optim_paras_info

    ! This mock object is required as we cannot simply pass in '' as it turns out.
    CHARACTER(225)                  :: file_sim_mock

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

    ! We intent to monitor the execution time of every evaluation of the criterion function.
    start = get_wtime()

    ! We allow for early termination due to maximum number of iterations or user request.
    IF (check_early_termination(maxfun, num_eval, crit_estimation)) THEN
        fort_criterion_scalar = HUGE_FLOAT
        RETURN
    END IF

    x_optim_free_unscaled = apply_scaling(x_optim_free_scaled, precond_matrix, 'undo')

    CALL construct_all_current_values(x_optim_all_unscaled, x_optim_free_unscaled, optim_paras, num_paras)

    CALL dist_optim_paras(optim_paras, x_optim_all_unscaled, dist_optim_paras_info)

    CALL fort_calculate_rewards_systematic(periods_rewards_systematic, num_periods, states_number_period, states_all, max_states_period, optim_paras)

    CALL fort_backward_induction(periods_emax, num_periods, is_myopic, max_states_period, periods_draws_emax, num_draws_emax, states_number_period, periods_rewards_systematic, mapping_state_idx, states_all, is_debug, is_interpolated, num_points_interp, edu_spec, optim_paras, file_sim_mock, .False.)

    CALL fort_contributions(contribs, periods_rewards_systematic, mapping_state_idx, periods_emax, states_all, data_est, periods_draws_prob, tau, num_periods, num_draws_prob, num_agents_est, num_obs_agent, num_types, edu_spec, optim_paras)

    fort_criterion_scalar = get_log_likl(contribs)

    IF (crit_estimation .OR. (maxfun == zero_int)) THEN

        num_eval = num_eval + 1

        CALL record_estimation(x_optim_free_scaled, x_optim_all_unscaled, fort_criterion_scalar, num_eval, num_paras, num_types, optim_paras, start)

        IF (dist_optim_paras_info .NE. zero_int) CALL record_warning(4)

    END IF

END FUNCTION
!******************************************************************************
!******************************************************************************
FUNCTION fort_criterion_parallel(x)

    !/* external objects    */

    REAL(our_dble), INTENT(IN)      :: x(:)
    REAL(our_dble)                  :: fort_criterion_parallel

    !/* internal objects    */

    REAL(our_dble)                  :: x_optim_all_unscaled(num_paras)
    REAL(our_dble)                  :: x_optim_free_unscaled(num_free)
    REAL(our_dble)                  :: contribs(num_agents_est)
    REAL(our_dble)                  :: start

    INTEGER(our_int), SAVE, ALLOCATABLE   :: num_states_slaves(:, :)
    INTEGER(our_int), SAVE, ALLOCATABLE   :: num_agents_slaves(:)
    INTEGER(our_int), SAVE, ALLOCATABLE   :: displs(:)

    INTEGER(our_int)                :: dist_optim_paras_info
    INTEGER(our_int)                :: i

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

#if MPI_AVAILABLE

    ! We intent to monitor the execution time of every evaluation of the criterion function.
    start = get_wtime()

    ! We allow for early termination due to maximum number of iterations or user request.
    IF (check_early_termination(maxfun, num_eval, crit_estimation)) THEN
        fort_criterion_parallel = HUGE_FLOAT
        RETURN
    END IF

    x_optim_free_unscaled = apply_scaling(x, precond_matrix, 'undo')

    CALL construct_all_current_values(x_optim_all_unscaled, x_optim_free_unscaled, optim_paras, num_paras)

    CALL MPI_Bcast(3, 1, MPI_INT, MPI_ROOT, SLAVECOMM, ierr)

    CALL MPI_Bcast(x_optim_all_unscaled, num_paras, MPI_DOUBLE, MPI_ROOT, SLAVECOMM, ierr)

    ! This extra work is only required to align the logging across the scalar and parallel implementation. In the case of an otherwise zero variance, we stabilize the algorithm. However, we want this indicated as a warning in the log file.
    CALL dist_optim_paras(optim_paras, x_optim_all_unscaled, dist_optim_paras_info)

    ! We need to know how the workload is distributed across the slaves.
    IF (.NOT. ALLOCATED(num_states_slaves)) THEN
        CALL distribute_workload(num_states_slaves, num_agents_slaves)

        ALLOCATE(displs(num_slaves))
        DO i = 1, num_slaves
            displs(i) = SUM(num_agents_slaves(:i - 1))
        END DO

    END IF

    contribs = -HUGE_FLOAT

    CALL MPI_GATHERV(contribs, 0, MPI_DOUBLE, contribs, num_agents_slaves, displs, MPI_DOUBLE, MPI_ROOT, SLAVECOMM, ierr)

    fort_criterion_parallel = get_log_likl(contribs)

    IF (crit_estimation .OR. (maxfun == zero_int)) THEN

        num_eval = num_eval + 1

        CALL record_estimation(x, x_optim_all_unscaled, fort_criterion_parallel, num_eval, num_paras, num_types, optim_paras, start)

        IF (dist_optim_paras_info .NE. zero_int) CALL record_warning(4)

    END IF

#endif

END FUNCTION
!******************************************************************************
!******************************************************************************
FUNCTION fort_dcriterion(x_optim_free_scaled)

    !/* external objects        */

    REAL(our_dble), INTENT(IN)      :: x_optim_free_scaled(:)
    REAL(our_dble)                  :: fort_dcriterion(SIZE(x_optim_free_scaled))

    !/* internals objects       */

    REAL(our_dble)                  :: ei(num_free)
    REAL(our_dble)                  :: d(num_free)
    REAL(our_dble)                  :: f0
    REAL(our_dble)                  :: f1

    INTEGER(our_int)                :: j

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

    ! Initialize containers
    ei = zero_dble

    ! Evaluate baseline
    f0 = criterion_function(x_optim_free_scaled)

    DO j = 1, num_free

        ei(j) = one_dble

        d = dfunc_eps * ei

        f1 = criterion_function(x_optim_free_scaled + d)

        fort_dcriterion(j) = (f1 - f0) / d(j)

        ei(j) = zero_dble

    END DO

END FUNCTION
!******************************************************************************
!******************************************************************************
SUBROUTINE construct_all_current_values(x_optim_all_unscaled, x_optim_free_unscaled, optim_paras, num_paras)

    !/* external objects        */

    INTEGER(our_int), INTENT(IN)        :: num_paras

    REAL(our_dble), INTENT(OUT)         :: x_optim_all_unscaled(num_paras)

    TYPE(OPTIMPARAS_DICT), INTENT(IN)   :: optim_paras

    REAL(our_dble), INTENT(IN)          :: x_optim_free_unscaled(COUNT(.not. optim_paras%paras_fixed))

    !/* internal objects        */

    INTEGER(our_int)                :: i
    INTEGER(our_int)                :: j

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

    j = 1

    DO i = 1, num_paras

        IF(optim_paras%paras_fixed(i)) THEN
            x_optim_all_unscaled(i) = x_all_start(i)
        ELSE
            x_optim_all_unscaled(i) = x_optim_free_unscaled(j)
            j = j + 1
        END IF

    END DO

END SUBROUTINE
!******************************************************************************
!******************************************************************************
SUBROUTINE get_precondition_matrix(precond_matrix, precond_spec, maxfun, x_optim_free_unscaled_start)

    !/* external objects    */

    REAL(our_dble), ALLOCATABLE, INTENT(OUT)    :: precond_matrix(:, :)

    TYPE(PRECOND_DICT), INTENT(IN)  :: precond_spec

    INTEGER(our_int), INTENT(IN)    :: maxfun

    REAL(our_dble), INTENT(IN)      :: x_optim_free_unscaled_start(num_free)

    !/* internal objects    */

    INTEGER(our_int)                :: i
    INTEGER(our_int)                :: u

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

    crit_estimation = .False.

    ALLOCATE(precond_matrix(num_free, num_free))

    CALL record_estimation(precond_matrix, x_optim_free_unscaled_start, optim_paras, .True.)

    IF ((precond_spec%type == 'identity') .OR. (maxfun == zero_int)) THEN
        precond_matrix = create_identity(num_free)
    ELSEIF (precond_spec%type == 'magnitudes') THEN
        precond_matrix = get_scales_magnitudes(x_optim_free_unscaled_start)
    ELSEIF (precond_spec%type == 'gradient') THEN
        precond_matrix = get_scales_gradient(x_optim_free_unscaled_start, precond_spec)
    ELSE
        STOP ' Not implemented ...'
    END IF

    ! Write out scaling matrix to allow for restart.
    5000 FORMAT(100000(1x,f45.15))

    OPEN(NEWUNIT=u, FILE='scaling.respy.out', ACTION='WRITE')

    DO i = 1, num_free
        WRITE(u, 5000) precond_matrix(i, :)
    END DO

    CLOSE(u)

    crit_estimation = .False.

END SUBROUTINE
!******************************************************************************
!******************************************************************************
FUNCTION get_scales_gradient(x_optim_free_start, precond_spec) RESULT(precond_matrix)

    !/* external objects    */

    REAL(our_dble)                              :: precond_matrix(num_free, num_free)

    REAL(our_dble), INTENT(IN)                  :: x_optim_free_start(num_free)

    TYPE(PRECOND_DICT), INTENT(IN)              :: precond_spec

    !/* internal objects    */

    REAL(our_dble)                              :: grad(num_free)
    REAL(our_dble)                              :: val

    INTEGER(our_int)                            :: i

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

    dfunc_eps = precond_spec%eps

    ! The precoditioning matrix needs to be initialized to approximate the gradient with a valid preconditioning matrix.
    precond_matrix = create_identity(num_free)
    grad = fort_dcriterion(x_optim_free_start)

    dfunc_eps = -HUGE_FLOAT

    precond_matrix = zero_dble

    DO i = 1, num_free

        val = ABS(grad(i))

        IF (val .LT. precond_spec%minimum) val = precond_spec%minimum

        precond_matrix(i, i) = val

    END DO

END FUNCTION
!******************************************************************************
!******************************************************************************
SUBROUTINE auto_adjustment_optimizers(optimizer_options, optimizer_used)

    !/* external objects    */

    TYPE(OPTIMIZER_COLLECTION), INTENT(INOUT)   :: optimizer_options

    CHARACTER(225), INTENT(IN)                  :: optimizer_used

    !/* internal objects    */

    INTEGER(our_int)                            :: npt

    REAL(our_dble)                              :: tmp(num_free)
    REAL(our_dble)                              :: rhobeg

    LOGICAL                                     :: is_misspecified

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

    IF (optimizer_used == 'FORT-NEWUOA') THEN

        npt = optimizer_options%newuoa%npt
        is_misspecified = (npt .LT. num_free + 2 .OR. npt .GT. ((num_free + 2) * (num_free + 1)) / 2)
        IF (is_misspecified) THEN
            optimizer_options%newuoa%npt = (2 * num_free) + 1
            CALL record_estimation('NEWUOA', optimizer_options%newuoa%npt)
        END IF

    END IF

    IF (optimizer_used == 'FORT-BOBYQA') THEN

        npt = optimizer_options%bobyqa%npt
        is_misspecified = (npt .LT. num_free + 2 .OR. npt .GT. ((num_free + 2) * (num_free + 1)) / 2)
        IF (is_misspecified) THEN
            optimizer_options%bobyqa%npt =  (2 * num_free) + 1
            CALL record_estimation('BOBYQA', optimizer_options%bobyqa%npt)
        END IF

        rhobeg = optimizer_options%bobyqa%rhobeg
        tmp = x_optim_bounds_free_scaled(2, :) - x_optim_bounds_free_scaled(1, :)

        rhobeg = optimizer_options%bobyqa%rhobeg
        is_misspecified = ANY(tmp .LT. rhobeg + rhobeg)
        IF (is_misspecified) THEN
            optimizer_options%bobyqa%rhobeg = MINVAL(tmp) * 0.5_our_dble
            optimizer_options%bobyqa%rhoend = optimizer_options%bobyqa%rhobeg * 1e-6
            CALL record_estimation('BOBYQA', optimizer_options%bobyqa%rhobeg, optimizer_options%bobyqa%rhoend)
        END IF

    END IF

END SUBROUTINE
!******************************************************************************
!******************************************************************************
END MODULE


! This file is part of fortran_tester
! Copyright 2015 Pierre de Buyl
! License: BSD

program test_tester_2
  use tester
  implicit none

  type(tester_t) :: test

  call test%init()

  call test%assert_equal(1, 2)

  call test%print()

end program test_tester_2


    program test_get_value_ptr

    use bmif_2_0, only: BMI_SUCCESS, BMI_FAILURE
    use bmiprmssoil
    use fixtures, only: config_file, status, print_1darray, print_i_1darray, &
        isreal4equalreal4, print_d_1darray

    implicit none

    type (bmi_prms_soil) :: m
    integer :: retcode

    retcode = test1()
    if (retcode.ne.BMI_SUCCESS) then
        stop BMI_FAILURE
    end if

    retcode = test2()
    if (retcode.ne.BMI_SUCCESS) then
        stop BMI_FAILURE
    end if

    contains

! Test getting r32 by nhru.
    function test1() result(code)
    character (len=*), parameter :: &
        var_name = "infil"
    integer, parameter :: rank = 1
    integer, parameter :: size = 14
    integer, parameter, dimension(rank) :: shape = (/ 14 /)
    real, parameter, dimension(shape(1)) :: &
        expected = (/0.02,0.02,0.03,0.04,0.02,0.01,0.03,0.01,0.04,0.09,0.05,0.08,0.05,0.06 /)
    real, pointer :: tptr(:)
    integer :: i, status, code

    double precision :: endtime

    status = m%initialize(config_file)
    status = m%get_end_time(endtime)
    do i = 1,int(endtime)
        status = m%update()
        if(i == endtime) then
            status = m%set_value(var_name, expected)
            status = m%get_value_ptr(var_name, tptr)
        endif
    enddo
    !status = m%get_value(var_name, tval)

    ! Visual inspection.
    write(*,*) "Test 1"
    write(*,*) "Expected"
    call print_1darray(expected, shape)
    write(*,*) "Get Value Ptr"
    call print_1darray(tptr, shape)

    code = BMI_SUCCESS
    do i = 1, shape(1)
        if (tptr(i).ne.expected(i)) then
            code = BMI_FAILURE
            exit
        end if
    end do
    status = m%finalize()

    end function test1
! Test getting r64 hru_area.
  function test2() result(code)
    character (len=*), parameter :: &
         var_name = "pfr_dunnian_flow"
    integer, parameter :: rank = 1
    integer, parameter :: size = 14
    integer, parameter, dimension(rank) :: shape = (/ 14 /)
    double precision, parameter, dimension(shape(1)) :: &
         expected = (/ 0.00D+000, 0.00D+000, 0.00D+000, 0.00D+000, &
                      0.00D+000, 0.00D+000, 0.00D+000, 0.00D+000, &
                      0.00D+000, 0.00D+000, 0.00D+000, 0.00D+000, &
                      0.00D+000, 0.00D+000 /)
    double precision, pointer :: tptr(:)

    integer :: i, status, code
    double precision :: endtime

    status = m%initialize(config_file)
    status = m%get_end_time(endtime)
    do i = 1,int(endtime)
        status = m%update()
        if(i == endtime) then
            status = m%get_value_ptr(var_name, tptr)
        endif
    enddo
    !status = m%get_value(var_name, tval)

    ! Visual inspection.
    write(*,*) "Test 2"
    write(*,*) "Expected"
    call print_d_1darray(expected, shape)
    write(*,*) "Get Value Ptr"
    call print_d_1darray(tptr, shape)

    code = BMI_SUCCESS
    do i = 1, size
       if (tptr(i).ne.expected(i)) then
          code = BMI_FAILURE
       end if
    end do
    status = m%finalize()
  end function test2

end program test_get_value_ptr


text:
!==============================================================================
!
!> @anchor  test_mean2true
!!
!> @brief   Program for testing transformation mean anomaly to true anomaly and vice versa.
!!		S/R from slam_astro_conversions are used
!!
!> @author  Eduard Gamper (EG)
!!
!> @date

!!
05.01.2015 (initial design)

!!

!!
!! @details Program for testing transformation mean anomaly to true anomaly and vice versa.
!!		Used equations in functions: Vallado, S.54 eq.2-6, S.56 eq. 2-13, 2-14
!!
!!
!------------------------------------------------------------------------

Program test_mean2true

use slam_astro_conversions
use slam_io, only: slam_message
use slam_types
use slam_time
use slam_orbit_types
use slam_math, only: eps15, eps9, eps6, eps3

implicit none

type(kepler_t) :: state_kepler			! State vector with Kepler elements
type(kepler_t) :: state_kepler_literature	! State vector with Kepler elements
						!              (1) = semimajor axis [km]
						!              (2) = eccentricity
						!              (3) = orbit inclination [rad]
						!              (4) = r.a. of asc. node [rad]
						!              (5) = argument of perigee [rad]
						!              (6) = mean anomaly [rad]

logical 	:: error = .false.		! Error flag
logical 	:: error_all = .false.		! Error flag

! QUELLE: http://www.tle.info/data/ISS_DATA.TXT
!
!               M50 Cartesian                         M50 Keplerian
!     -----------------------------------       --------------------------------
!      X    =        3831218.94                 A    =         6778927.15  meter
!      Y    =       -3769807.71  meter          E    =           .0013488
!      Z    =        4131046.64                 I    =           51.67403
!      XDOT =       6287.659929                 Wp   =           38.31673
!      YDOT =       2227.041971  meter/sec      RA   =          172.87570  deg
!      ZDOT =      -3782.037544                 TA   =           90.71589
!                                               MA   =           90.56134
!                                               Ha   =            220.633  n.mi
!                                               Hp   =            212.869

! subroutine mean2true(         &
!                       Ecc,    &  ! <-- DBL  eccentricity
!                       M,      &  ! <-- DBL  mean anomaly / rad
!                       E0,     &  ! --> DBL  eccentric anomaly / rad
!                       Nu      &  ! --> DBL  true anomaly / rad
!                     )

! subroutine true2mean(       &
!                       Ecc,  & ! <-- DBL eccentricity
!                       Nu,   & ! <-- DBL true anomly in rad
!                       E0,   & ! --> DBL eccentric anomaly in rad
!                       M     & ! --> DBL mean anomaly in rad
!                     )

state_kepler_literature%ecc  = 0.0013488d0
state_kepler_literature%man  = 90.56134d0/180*pi
state_kepler_literature%tran = 90.71589d0/180*pi

state_kepler%ecc = state_kepler_literature%ecc
state_kepler%man = state_kepler_literature%man

write(*,*) ''
write(*,*) 'Testing transformation mean anomaly to true anomaly and vice versa'
write(*,*) 'Functions mean2true() and true2mean() from module slam_astro_legacy are used'
write(*,*) ''
write(*,*) '------- Example ISS -------'
write(*,*) '--Data from literature (ISS-Orbit):'
write(*,*) 'Eccentricity: ', state_kepler_literature%ecc
write(*,*) 'Mean anomaly: ', state_kepler_literature%man*180/pi
write(*,*) 'True anomaly: ', state_kepler_literature%tran*180/pi

!Transform mean anomaly to true anomaly
CALL mean2true(state_kepler%ecc, state_kepler%man, state_kepler%ecan, state_kepler%tran)

write(*,*) '--Transform mean anomaly to true anomaly (ISS-Orbit):'
write(*,*) 'True anomaly: ', state_kepler%tran*180/pi
write(*,*) 'Eccentric anomaly: ', state_kepler%ecan*180/pi

!** Compare results of s/r mean2true(); use of eps6 because iteration from mean 2 true has allowed deviation of eps9
if (abs(state_kepler%tran - state_kepler_literature%tran) .gt. eps6) error = .true.

!** Display conclusion
if (error) then
  CALL slam_message('mean2true not passed.',1)
  error_all = .true.
else
  CALL slam_message('mean2true passed.', 1)
end if

!Transform true anomaly back to mean anomaly
CALL true2mean(state_kepler%ecc, state_kepler%tran, state_kepler%ecan, state_kepler%man)

write(*,*) '--Transform true anomaly back to mean anomaly (ISS-Orbit):'
write(*,*) 'Mean anomaly: ', state_kepler%man*180/pi
write(*,*) 'Eccentric anomaly: ', state_kepler%ecan*180/pi

!** Set error flag back to false
error = .false.

!** Compare results of s/r true2mean()
if (abs(state_kepler%man - state_kepler_literature%man) .gt. eps6) error = .true.

!** Display conclusion
if (error) then
  CALL slam_message('true2mean not passed.',1)
  error_all = .true.
else
  CALL slam_message('true2mean passed.', 1)
end if

! Example: High eccentricity
! 0 MOLNIYA 3-50
! 1 25847U 99036A   15310.12862371 -.00000116 +00000-0 -13949-3 0  9992
! 2 25847 062.6934 338.4194 7400458 260.6635 016.9238 02.00678845119668

! Mean Anomaly: 016.9238 deg
! Eccentricity: 0.7400458

state_kepler_literature%ecc   = 0.7400458d0
state_kepler_literature%man   = 16.9238d0/180*pi
state_kepler_literature%tran  = 99.203171087928723d0/180*pi !This is a calculated value (by s/r mean2true)

state_kepler%ecc = state_kepler_literature%ecc
state_kepler%man = state_kepler_literature%man
! state_kepler%tran=90.71589d0/180*pi

write(*,*) ''
write(*,*) '------- Example Molniya -------'
write(*,*) '--Data from literature (Molniya-Orbit):'
write(*,*) 'Eccentricity: ', state_kepler%ecc
write(*,*) 'Mean anomaly: ', state_kepler%man*180/pi
! write(*,*) 'True anomaly: ', state_kepler%tran*180/pi

!Transform mean anomaly to true anomaly (Molniya-Orbit)
CALL mean2true(state_kepler%ecc, state_kepler%man, state_kepler%ecan, state_kepler%tran)

write(*,*) '--Transform mean anomaly to true anomaly (Molniya-Orbit):'
write(*,*) 'True anomaly: ', state_kepler%tran*180/pi
write(*,*) 'Eccentric anomaly: ', state_kepler%ecan*180/pi

!** Set error flag back to false
error = .false.

!** Compare results of s/r mean2true()
if (abs(state_kepler%tran - state_kepler_literature%tran) .gt. eps6) error = .true.
write(*,*) abs(state_kepler%tran - state_kepler_literature%tran)

!** Display conclusion
if (error) then
  CALL slam_message('mean2true not passed.',1)
  error_all = .true.
else
  CALL slam_message('mean2true passed.', 1)
end if

!Transform true anomaly back to mean anomaly (Molniya-Orbit)
CALL true2mean(state_kepler%ecc, state_kepler%tran, state_kepler%ecan, state_kepler%man)

write(*,*) '--Transform true anomaly back to mean anomaly (Molniya-Orbit):'
write(*,*) 'Mean anomaly: ', state_kepler%man*180/pi
write(*,*) 'Eccentric anomaly: ', state_kepler%ecan*180/pi

!** Set error flag back to false
error = .false.

!** Compare results of s/r true2mean()
if (abs(state_kepler%man - state_kepler_literature%man) .gt. eps6) error = .true.

!** Display conclusion
if (error) then
  CALL slam_message('true2mean not passed.',1)
  error_all = .true.
else
  CALL slam_message('true2mean passed.', 1)
end if

!** Display conclusion for complete test
if (error_all) then
  CALL slam_message('Failed.',1)
  error_all = .true.
else
  CALL slam_message('All tests passed.', 1)
end if

end Program test_mean2true

max_stars_repo_name: Space-Systems/libslam

max_issues_repo_name: Space-Systems/libslam

max_forks_repo_name: Space-Systems/libslam


text:
      DOUBLE COMPLEX FUNCTION ZDOTC( N, X, INCX, Y, INCY )

      INTEGER            INCX, INCY, N
      DOUBLE COMPLEX     X( * ), Y( * )

      DOUBLE COMPLEX     ZERO
      PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ) )

      DOUBLE COMPLEX     ZTEMP
      EXTERNAL           ZDOTC_SUB

      IF( N.GT.0 ) THEN
         CALL ZDOTC_SUB( N, X, INCX, Y, INCY, ZTEMP )
         ZDOTC = ZTEMP
      ELSE
         ZDOTC = ZERO
      END IF

      RETURN
      END

ext: f

max_stars_repo_name: stip/FLENS

max_issues_repo_name: stip/FLENS

max_forks_repo_name: stip/FLENS