text: #redirect New Pariss Nails

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!
! m_exception --
!
!   Provides services to generate different levels of exceptions
!   and display the message of the exception.
!   Five levels of exceptions can be generated, from the lowest
!   to the highest level :
!   - information, warning : just print a message and continue
!   - error, fatal_error and failure : prints a message and stop the execution
!
! OVERVIEW
!
!   Simple use-case
!
!   Suppose that one would like to compute the square root of one
!   real value. The "compute_sqrt" function takes one positive real argument,
!   and if the argument is negative, one cannot compute the square root so
!   that one would generate an error. In the following example, extracted from the
!   unit tests included in the project, one uses the static method "exception_raiseError"
!   to display a user-friendly message and stop the execution of the program
!
!    function compute_sqrt ( value ) result ( root )
!      use m_exception
!      implicit none
!      real, intent(in) :: value
!      real :: root
!      if ( value < 0. ) then
!        call exception_raiseError ( "Value is negative in compute_sqrt" )
!      else
!        root = sqrt ( value )
!      endif
!   end function compute_sqrt
!
!   real :: root
!   root = compute_sqrt ( -1. )
!
!   In the previous example, the standard output is written so that the
!   following message appears on screen :
!
!   Error.
!   Message: Value is negative in compute_sqrt
!
!   Controlling the execution
!
!   The client code can control the behaviour of the component each time
!   an exception is raised.
!   The default behaviour is to stop the execution. This can be modified
!   by calling "exception_setstoponerror" in order to continue the execution,
!   even if error, fatal error or failure exceptions are raised.
!
!   In the following example, the static method "exception_setstoponerror" is
!   called so that an error does not interrupt the execution.
!
!   call exception_setstoponerror ( .false. )
!   call exception_raiseError ( "There is an error, but the execution will continue." )
!
! Controlling output
!
!   The default behaviour is to write messages onto the standard output
!   each time an exception is raised.
!   This can be modified in two ways :
!   * the first possibility is to disable the writing of the messages
!     with "exception_logactive". This feature might be useful in the case
!     where a component has known bugs but generates lots of unwanted
!     exceptions messages.
!   * the second possibility is to connect the component to an existing unit
!     with "exception_setlogunit", so that the messages are written
!     on the given logical unit number.
!     This allows for example to write on an existing log file, may be the log
!     file manage by the m_logger component included in the project.
!
!   In the following example, the client code first disables all output,
!   set "stoponerror" to false and generates an error which is not displayed
!   and does not interrupt the execution.
!
!   call exception_setstoponerror ( .false. )
!   call exception_logactive ( .false. )
!   call exception_raiseError ( "This message will not be displayed and the execution will continue." )
!   call exception_logactive ( .true. )
!   call exception_raiseError ( "This message WILL be displayed and the execution will continue." )
!
!   In the following example, the client code connects the m_exception component to
!   an existing unit so that the exception messages are written onto a client log file.
!
!     log_fileunit = 12
!     call exception_setstoponerror ( .false. )
!     open ( log_fileunit , FILE= "log_file.log" )
!     call exception_setlogunit ( log_fileunit )
!     call exception_raiseError ( "This message will be written in log_file.log and the execution will continue." )
!     call exception_setlogunit ( 0 )
!     call exception_raiseError ( "This message will be written on standard output and the execution will continue." )
!     close ( log_fileunit )
!
!   In the following example, the client code connects the m_exception component to
!   the logfile manage by m_logger. This way, the exception messages are collected in the
!   unique log file of the client code.
!
!     call log_startup ( "log_file.log" , append=.true. )
!     call log_cget ( "logfileunit" , log_fileunit )
!     call exception_setstoponerror ( .false. )
!     call exception_setlogunit ( log_fileunit )
!     call exception_raiseError ( "This message will be written in log_file.log and the execution will continue." )
!     call log_shutdown ()
!
! Pseudo-catch
!
!   The client code can use a pseudo-catch system which provides
!   a simple way to manage exceptions which are raised at a lower
!   level in the call stack. This allows to provide special
!   treatments when exceptions are generated, without modifiying
!   all lower level subroutines/function, but just by inserting
!   exception management when needed.
!   Suppose that you have a subroutine which source code is :
!     subroutine yoursubroutine ()
!       use m_exception, only : exception_raiseFatalError
!       implicit none
!       [...]
!       call exception_raiseFatalError ( "Wrong blabla in yoursubroutine" )
!       [...]
!     end subroutine yoursubroutine
!   When calling the subroutine "yoursubroutine", one may wonder if exceptions
!   have been generated so that these errors may be processed, or not.
!   One can use the "exception_catch" service to compute the status
!   of one subroutine and manage that status :
!     use m_exception, only : exception_catch, &
!         EXCEPTION_INFORMATION, &
!         EXCEPTION_WARNING &
!         EXCEPTION_ERROR &
!         EXCEPTION_FATAL_ERROR &
!         EXCEPTION_FAILURE
!     integer :: status
!     call exception_catch ( yoursubroutine , status )
!     select case ( status )
!     case ( EXCEPTION_INFORMATION )
!        write(6,*) "Information"
!     case ( EXCEPTION_WARNING )
!        write(6,*) "Warning"
!     case ( EXCEPTION_ERROR , EXCEPTION_FATAL_ERROR , EXCEPTION_FAILURE )
!        write(6,*) "Fatal error"
!     case default
!        write(6,*) "No problem, continue."
!     end select
!
! TODO
!   - design a more powerful exception management system, which manages exceptions
!   through the call stack
!
! Copyright (c) 2008 Michael Baudin
!
! $Id: m_exception.f90,v 1.5 2013/05/13 08:03:21 knystrom Exp $
!
module m_exception
  implicit none
  private
  !
  ! Public methods
  !
  public :: exception_getcounter
  public :: exception_getlogunit
  public :: exception_initcounter
  public :: exception_islogactive
  public :: exception_logactive
  public :: exception_raiseError
  public :: exception_raiseFailure
  public :: exception_raiseFatalError
  public :: exception_raiseInformation
  public :: exception_raiseWarning
  public :: exception_report
  public :: exception_setlogunit
  public :: exception_setstoponerror
  public :: exception_catch
  !
  ! Tags to manage error
  !
  integer, parameter, public :: EXCEPTION_OK = 0
  integer, parameter, public :: EXCEPTION_INFORMATION = 1
  integer, parameter, public :: EXCEPTION_WARNING = 2
  integer, parameter, public :: EXCEPTION_ERROR = 3
  integer, parameter, public :: EXCEPTION_FATAL_ERROR = 4
  integer, parameter, public :: EXCEPTION_FAILURE = 5
  integer, parameter, public :: EXCEPTION_SIZE = 5
  !
  ! Set to true to stop the execution after a stoping error.
  !
  logical :: exception_stoponerror = .true.
  !
  ! Counters for the different levels of errors
  !
  integer :: exception_nbinformation = 0
  integer :: exception_nbwarning = 0
  integer :: exception_nberror = 0
  integer :: exception_nbfatalerror = 0
  integer :: exception_nbfailure = 0
  !
  ! Maximum number of characters in a message generated by the exception
  !
  integer , parameter :: EXCEPTION_MAXIMUM_LENGTH = 500
  !
  ! Set exception_log_active to true to activate the exception logging
  ! to the unit # exception_log_unit
  !
  integer :: exception_log_unit = 0
  logical :: exception_log_unit_active = .false.
  !
  ! Set exception_log_active to true to display the exception messages
  !
  logical :: exception_log_active = .true.
contains
  !
  ! exception_raiseFatalError --
  !   Generates a fatal error message and stops the execution.
  !
  subroutine exception_raiseFatalError ( message )
    character(len=*), intent(in) :: message
    exception_nbfatalerror = exception_nbfatalerror + 1
    call exception_print(EXCEPTION_FATAL_ERROR, message )
    call exception_stop ()
  end subroutine exception_raiseFatalError
  !
  ! exception_raiseError --
  !   Generates an error message and stops the execution
  !
  subroutine exception_raiseError ( message )
    character(len=*), intent(in) :: message
    exception_nberror = exception_nberror + 1
    call exception_print(EXCEPTION_ERROR, message )
    call exception_stop ()
  end subroutine exception_raiseError
  !
  ! exception_raiseWarning --
  !   Generates a warning message
  !
  subroutine exception_raiseWarning ( message )
    character(len=*), intent(in) :: message
    exception_nbwarning = exception_nbwarning + 1
    call exception_print(EXCEPTION_WARNING, message )
  end subroutine exception_raiseWarning
  !
  ! exception_raiseInformation --
  !   Generates an information message.
  !
  subroutine exception_raiseInformation ( message )
    character(len=*), intent(in) :: message
    exception_nbinformation = exception_nbinformation + 1
    call exception_print(EXCEPTION_INFORMATION, message )
  end subroutine exception_raiseInformation
  !
  ! exception_raiseFailure --
  !   Generates a failure message and stops the execution
  !
  subroutine exception_raiseFailure ( message )
    character(len=*), intent(in) :: message
    exception_nbfailure = exception_nbfailure + 1
    call exception_print(EXCEPTION_FAILURE, message )
    call exception_stop ()
  end subroutine exception_raiseFailure
  !
  ! exception_print --
  !   Prints a message on the standard output about
  !   the current exception type, the origin of the error and the
  !   message of the error.
  ! Arguments :
  !   exceptionCode : an integer exception code between : EXCEPTION_INFORMATION,
  !     EXCEPTION_WARNING, EXCEPTION_ERROR, EXCEPTION_FATAL_ERROR, EXCEPTION_FAILURE
  !   origin : the name of the element generating the exception
  !   message : a message to display
  !
  subroutine exception_print ( exceptionCode , message )
    integer, intent(in) :: exceptionCode
    character(len=*), intent(in) :: message
    character ( len = EXCEPTION_MAXIMUM_LENGTH ) :: msg
    select case (exceptionCode)
    case (EXCEPTION_INFORMATION)
       write ( msg , * ) "Information."
    case (EXCEPTION_WARNING)
       write(msg,*) "Warning."
    case (EXCEPTION_ERROR)
       write(msg,*) "Error."
    case (EXCEPTION_FATAL_ERROR)
       write(msg,*) "Fatal error."
    case (EXCEPTION_FAILURE)
       write(msg,*) "Failure."
    case default
       write(msg,*) "Bad value for the exception code."
    end select
    call exception_logmsg ( msg )
    write(msg,*) "Message: ", trim(message)
    call exception_logmsg ( msg )
  end subroutine exception_print
  !
  ! exception_stop --
  !   Stop the execution after an exception has been raised.
  ! Note :
  !   In debug mode, put a breakpoint here allows to see the exception
  !   before stopping.
  !
  subroutine exception_stop ()
    if ( exception_stoponerror ) then
       STOP
    endif
  end subroutine exception_stop
  !
  ! exception_setstoponerror --
  !   Configure the component to manage the behaviour when an exception
  !   is generated.
  ! Arguments :
  !   stoponerror : if true, then when an error, a fatal_error or a failure
  !     is generated, then stop the execution.
  !     If .false., then the same kinds of exceptions do not stop the
  !     execution
  !
  subroutine exception_setstoponerror ( stoponerror )
    logical, intent(in) :: stoponerror
    exception_stoponerror = stoponerror
  end subroutine exception_setstoponerror
  !
  ! exception_logmsg --
  !   Write the given message onto the default unit.
  !
  subroutine exception_logmsg ( message )
    character(len=*), intent(in) :: message
    if ( exception_log_active ) then
       if ( exception_log_unit_active ) then
          write(exception_log_unit,*) trim(message)
       else
          write(*,*) trim(message)
       endif
    endif
  end subroutine exception_logmsg
  !
  ! exception_initcounter --
  !   Reset to 0 all the exceptions counters
  !
  subroutine exception_initcounter ()
    exception_nbinformation = 0
    exception_nbwarning = 0
    exception_nberror = 0
    exception_nbfatalerror = 0
    exception_nbfailure = 0
  end subroutine exception_initcounter
  !
  ! exception_getcounter --
  !   Returns an array of size EXCEPTION_SIZE indicating the number
  !   of type of exceptions.
  ! Arguments :
  !   excepcounters ( 1:EXCEPTION_SIZE) : array of integers
  !     excepcounters ( iexcept ) is the total number of exceptions
  !     of type #iexcept generated since the begining of the execution.
  !     The possible values for iexcept are the following.
  !     - counter ( EXCEPTION_INFORMATION ) : total number of information exceptions raised
  !     - counter ( EXCEPTION_WARNING ) : total number of warning exceptions raised
  !     - counter ( EXCEPTION_ERROR ) : total number of error exceptions raised
  !     - counter ( EXCEPTION_FATAL_ERROR ) : total number of fatal error exceptions raised
  !     - counter ( EXCEPTION_FAILURE ) : total number of failure exceptions raised
  !
  subroutine exception_getcounter ( counter )
    integer, dimension(1:EXCEPTION_SIZE), intent(out) :: counter
    counter ( EXCEPTION_INFORMATION ) = exception_nbinformation
    counter ( EXCEPTION_WARNING ) = exception_nbwarning
    counter ( EXCEPTION_ERROR ) = exception_nberror
    counter ( EXCEPTION_FATAL_ERROR ) = exception_nbfatalerror
    counter ( EXCEPTION_FAILURE ) = exception_nbfailure
  end subroutine exception_getcounter
  !
  ! exception_report --
  !   Writes on the current exception unit number a report
  !   which details the number of exceptions of all types
  !   since the begining of the execution or the last reset of
  !   all counters.
  !
  subroutine exception_report ()
    character ( len = 300 ) :: message
    call exception_logmsg ( "Exception report" )
    write ( message , * ) "Number of informations :" , exception_nbinformation
    call exception_logmsg ( message )
    write ( message , * ) "Number of warnings :" , exception_nbwarning
    call exception_logmsg ( message )
    write ( message , * ) "Number of errors :" , exception_nberror
    call exception_logmsg ( message )
    write ( message , * ) "Number of fatal errors :" , exception_nbfatalerror
    call exception_logmsg ( message )
    write ( message , * ) "Number of failures :" , exception_nbfailure
    call exception_logmsg ( message )
  end subroutine exception_report
  !
  ! exception_setlogunit --
  !   Set the unit number onto which the messages are output.
  !   If the unitnumber is 0, the messages are written
  !   to the standard output (unit *).
  ! Arguments :
  !   unitnumber : the unit number
  !
  subroutine exception_setlogunit ( unitnumber )
    integer, intent(in) :: unitnumber
    if ( unitnumber == 0 ) then
       exception_log_unit_active = .false.
       exception_log_unit = unitnumber
    else
       exception_log_unit_active = .true.
       exception_log_unit = unitnumber
    endif
  end subroutine exception_setlogunit
  !
  ! exception_getlogunit --
  !   Returns the positive unit number onto which the messages are output,
  !   if enabled or a negative integer if the feature is disabled.
  ! Arguments :
  !
  !
  function exception_getlogunit ( ) result ( logunit )
    integer :: logunit
    logunit = exception_log_unit
  end function exception_getlogunit
  !
  ! exception_logactive --
  !   If the boolean argument bool is true, enable the logging
  !   of the exceptions messages.
  !   If the boolean argument bool is false, disable the logging
  !   of the exceptions messages.
  ! Arguments :
  !   bool : if true, active the logging of the messages
  !
  subroutine exception_logactive ( bool )
    logical, intent(in) :: bool
    exception_log_active = bool
  end subroutine exception_logactive
  !
  ! exception_islogactive --
  !   Returns .true. if the current exception messages are written,
  !   either on standard output or into a log file.
  !
  function exception_islogactive ( ) result ( islogactive )
    logical :: islogactive
    islogactive = exception_log_active
  end function exception_islogactive
  !
  ! exception_catch --
  !   Calls the given subroutine and returns the integer associated
  !   with last exception, higher level, code or 0 if no exception
  !   was raised.
  ! Argument :
  !   callback : the subroutine to call back
  !   status : the integer corresponding to the last exception, if any.
  ! Caution !
  !   The internal algorithm is based on the exception counters,
  !   which implies that any call to exception_initcounter in the client
  !   code can result in a wrong status.
  !
  subroutine exception_catch ( callback , status )
    interface interfacecallback
       subroutine callback ()
       end subroutine callback
    end interface interfacecallback
    integer, intent(out) :: status
    integer, dimension(1:EXCEPTION_SIZE) :: counter_before
    integer, dimension(1:EXCEPTION_SIZE) :: counter_after
    integer :: exception_code
    !
    ! Get the counters before
    !
    call exception_getcounter ( counter_before )
    !
    ! Call the callback
    !
    call callback ()
    !
    ! Get the counters after
    !
    call exception_getcounter ( counter_after )
    !
    ! Compare before / after and set the corresponding status.
    !
    ! Process the loop, from the highest level of exception to the lowest.
    status = EXCEPTION_OK
    do exception_code = EXCEPTION_SIZE , 1 , -1
       if ( counter_after ( exception_code ) > counter_before ( exception_code ) ) then
          status = exception_code
          exit
       endif
    enddo
  end subroutine exception_catch
end module m_exception

max_stars_repo_name: philip-peterson/flibs

max_issues_repo_name: philip-peterson/flibs

max_forks_repo_name: philip-peterson/flibs


text:
      Module mba3d_auxSF
C
      use mba3d_auxSE
      use mba3d_auxSR
      use mba3d_error
      use mba3d_findSE
      use mba3d_makQ
C
      contains
C
C ================================================================
      Subroutine makSF(
C ================================================================
     &           iP, lF, iFs, lE, iEs, IPEs, IPF, IFE, MaxS,
     &           lP1, iP1s, lP2, iP2s, ICP2s,
     &           lS, IPSs, IESs, lDF, iDFs, lDE, iDEs,
     &           flagFACE, flagEDGE)
C ================================================================
C Compute:
C   (a) list of faces IPSs(3, lS) around point iP is created.
C
C   (b) list of surface points connected with point iP
C
C   (c) edge case : list of points on an edge passing through iP
C       face case : list of points on the face passing through iP
C ================================================================
      Integer iFs(*), iEs(*)
      Integer IPEs(5, *), IPF(4, *), IFE(4, *)
      Integer iP1s(*), iP2s(*), ICP2s(*)
      Integer IPSs(3, *), IESs(*), iDFs(*), iDEs(*)
      Logical flagFACE, flagEDGE

C group (Local variables)
      Integer iref(5), iPt(3)
      Logical flag

C ================================================================
      iref(1) = 1
      iref(2) = 2
      iref(3) = 3
      iref(4) = 4
      iref(5) = 1

      lP1 = 0
      lP2 = 0

      lS  = 0

      lDF = 0
      lDE = 0

      flagFACE = .FALSE.
      flagEDGE = .FALSE.

      Do n = 1, lE
         flag = .FALSE.
         Do 10 i1 = 1, 4
            If(IPEs(i1, n).EQ.iP) Then
               flag = .TRUE.

               lDE = lDE + 1
               iDEs(lDE) = n

               i2 = i1
               Do i = 1, 3
                  i2 = iref(i2 + 1)
                  iPt(i) = IPEs(i2, n)

                  Call findSE(lP2, iP2s, iPt(i), nPt)

                  If(nPt.EQ.0) Then
                     lP2 = lP2 + 1
                     If(lP2.GT.MaxS) Goto 1000

                     iP2s(lP2) = iPt(i)
                     ICP2s(lP2) = 0
                  End if
               End do

               Do m = 1, lS
                  If(check33(iPt(1), iPt(2), iPt(3),
     &               IPSs(1, m), IPSs(2, m), IPSs(3, m))) Goto 10
               End do

               lS = lS + 1
               If(lS.GT.MaxS) Goto 1000

               Do i = 1, 3
                  IPSs(i, lS) = iPt(i)
               End do
               IESs(lS) = n

               Goto 15
            End if
 10     Continue

 15     If(flag) Then
           Do 20 i1 = 1, 4
              iFt = IFE(i1, iEs(n))

              If(iFt.NE.0) Then
                 i2 = iref(i1 + 1)
                 i3 = iref(i2 + 1)

                 iPt(1) = IPEs(i1, n)
                 iPt(2) = IPEs(i2, n)
                 iPt(3) = IPEs(i3, n)

                 If(check13(iP, iPt(1), iPt(2), iPt(3))) Then
                    iCLRs = IPF(4, iFt)
                    If(iCLRs.LT.0) Call errMes(6001, 'auxSF',
     &                                        'system error')

                    Do i = 1, 3
                       Call findSE(lP2, iP2s, iPt(i), nPt)
                       If(nPt.GT.0) Then
                          If(ICP2s(nPt).EQ.0) Then
                             flagFACE = .TRUE.

                             ICP2s(nPt) = iCLRs
                          Else If(ICP2s(nPt).NE.iCLRs) Then
                             flagEDGE = .TRUE.

                             Call findSE(lP1, iP1s, iPt(i), mPt)
                             If(mPt.EQ.0) Then
                                lP1 = lP1 + 1
                                If(lP1.GT.MaxS) Goto 1000

                                iP1s(lP1) = iPt(i)
                             End if

                             ICP2s(nPt) = ICP2s(nPt) + iCLRs
                           End if
                       End if
                    End do

                    Call findSE(lF, iFs, iFt, nFt)
                    Call findSE(lDF, iDFs, nFt, k)
                    If(k.LE.0) Then
                       lDF = lDF + 1
                       If(lDF.GT.MaxS) Goto 1000

                       iDFs(lDF) = nFt
                    End if
                 End if
              End if
 20        Continue
        End if
      End do

      If(flagEDGE) Then
         Goto 9000
      Else If(flagFACE) Then
         Do n = 1, lP2
            If(ICP2s(n).NE.0) Then
               lP1 = lP1 + 1
               iP1s(lP1) = iP2s(n)
            End if
         End do
      Else
         lP1 = lP2
         Do n = 1, lP1
            iP1s(n) = iP2s(n)
         End do
      End if

 9000 Return
 1000 Call errMes(1007, 'makSF', 'local variable MaxS is small')
      End Subroutine makSF

C ================================================================
      Subroutine shutF(ZZ1, ZZ2, xyp1, xyp2, xyp3, flag)
C ================================================================
C Routine returns .TRUE. if the ray [ZZ1, ZZ2) intersects
C triangle {xyz1, xyz2, xyz3}. If flag = .TRUE. the triangle is
C assumed to be the open domain.
C ================================================================
      Real*8  ZZ1(3),  ZZ2(3)
      Real*8  xyp1(3), xyp2(3), xyp3(3)
      Logical flag

C group (Local variables)
      Real*8  v1, v2
      Logical flagOPEN

C ================================================================
      flagOPEN = flag
      flag = .FALSE.

      v1 = calVol(ZZ1, ZZ2, xyp3, xyp1)
      v2 = calVol(ZZ1, ZZ2, xyp3, xyp2)

      If(flagOPEN) Then
         If(v1 * v2.GE.0D0) Goto 1000
      Else
         If(v1 * v2.GT.0D0) Goto 1000
      End if

      v1 = calVol(ZZ1, ZZ2, xyp2, xyp1)
      v2 = calVol(ZZ1, ZZ2, xyp2, xyp3)

      If(flagOPEN) Then
         If(v1 * v2.GE.0D0) Goto 1000
      Else
         If(v1 * v2.GT.0D0) Goto 1000
      End if

      v1 = calVol(ZZ1, ZZ2, xyp1, xyp2)
      v2 = calVol(ZZ1, ZZ2, xyp1, xyp3)

      If(flagOPEN) Then
         If(v1 * v2.GE.0D0) Goto 1000
      Else
         If(v1 * v2.GT.0D0) Goto 1000
      End if

      flag = .TRUE.
 1000 Return
      End Subroutine shutF

C ================================================================
      Real*8 Function projectF(xyp, xy1, xy2, xy3, xyt)
C ================================================================
C Routine finds a point xyt lying inside triangular face
C {xy1, xy2, xy3} which is the closest point to xyp. It returns
C the distance to this point.
C
C *** Remarks:
C        1. We minimize the distance ||a * x21 + b * x31 - p||^2
C           w.r.t. a and b.
C ================================================================
      Real*8  xyp(3), xy1(3), xy2(3), xy3(3), xyt(3)

c group (Local variables)
      Real*8  v1(3), v2(3), vp(3)
      Real*8  m11, m12, m22, f1, f2, det, a, b

C ================================================================
      Do i = 1, 3
         v1(i) = xy2(i) - xy1(i)
         v2(i) = xy3(i) - xy1(i)
         vp(i) = xyp(i) - xy1(i)
      End do

      m11 = DotMul(v1, v1)
      m12 = DotMul(v1, v2)
      m22 = DotMul(v2, v2)

      f1 = DotMul(vp, v1)
      f2 = DotMul(vp, v2)

c ... solve a linear system (Wronskian is always non-zero)
      det = m22 * m11 - m12 * m12

      a = (m22 * f1 - m12 * f2) / det
      b = (m11 * f2 - m12 * f1) / det

c ... analyze few cases (seach along b, a and a + b = 1)
      If(a.LT.0D0) Then
         a = 0D0
         b = max(0D0, f2 / m22)
         b = min(1D0, b)
      Else If(b.LT.0D0) Then
         b = 0D0
         a = max(0D0, f1 / m11)
         a = min(1D0, a)
      Else If(a + b.GT.1D0) Then
         a = max(0D0, (f1 - m12) / (m11 - m12))
         a = min(1D0, a)
         b = 1D0 - a
      End if

c ... compute point the closest point
      Do i = 1, 3
         xyt(i) = xy1(i) + a * v1(i) + b * v2(i)
      End do

      projectF = calEdge(xyp, xyt)

      Return
      End Function projectF

C ================================================================
      Subroutine bariCoords(xyp, xy1, xy2, xy3, b1, b2, b3)
C ================================================================
C Routine computes barycentric coordinates of point xyp lying
C inside triangle {xy1, xy2, xy3}. They coordinates are always
C positive and sum to 1 upto round-off errors.
C ================================================================
      Real*8  xyp(3), xy1(3), xy2(3), xy3(3)
      Real*8  b1, b2, b3

c group (Local variables)
      Real*8  a(3), b(3), c(3), art, ar1

C ================================================================
c ... square of area of the triangle
      Do i = 1, 3
         a(i) = xy2(i) - xy1(i)
         b(i) = xy3(i) - xy1(i)
      End do

      Call VecMul(a, b, c)
      art = c(1) ** 2 + c(2) ** 2 + c(3) ** 2

c ... first baricentric coordinate
      Do i = 1, 3
         a(i) = xy2(i) - xyp(i)
         b(i) = xy3(i) - xyp(i)
      End do

      Call VecMul(a, b, c)
      ar1 = c(1) ** 2 + c(2) ** 2 + c(3) ** 2

      b1 = dsqrt(ar1 / art)

c ... second baricentric coordinate
      Do i = 1, 3
         a(i) = xy1(i) - xyp(i)
      End do

      Call VecMul(a, b, c)
      ar1 = c(1) ** 2 + c(2) ** 2 + c(3) ** 2

      b2 = dsqrt(ar1 / art)

c ... third baricentric coordinate
      b3 = 1 - b2 - b1
      If(b1 + b2.GT.1D0) Then
         b1 = min(1D0, b1)
         b2 = max(0D0, 1 - b1)
         b3 = 0D0
      End if

      Return
      End Subroutine bariCoords

C ================================================================
      Real*8 Function angle2Faces( XYA, XYB, XYC, XYD )
C ================================================================
C Cosine of angle between two faces {A,B,C} and {A,B,D}.
C
C Remark: the angle is equal to -1 for flat faces.
C ================================================================
      Real*8  XYA(3), XYB(3), XYC(3), XYD(3)
      Real*8  V1(3), V2(3), V3(3), XYN(3), XYM(3)
      Real*8  d1, d2

      Do i = 1, 3
         V1(i) = XYB(i) - XYA(i)
         V2(i) = XYC(i) - XYB(i)
         V3(i) = XYD(i) - XYB(i)
      End do

      Call VecMul(V1, V2, XYN)
      Call VecMul(V1, V3, XYM)

      d1 = XYN(1) ** 2 + XYN(2) ** 2 + XYN(3) ** 2
      d2 = XYM(1) ** 2 + XYM(2) ** 2 + XYM(3) ** 2

      angle2Faces = (XYN(1) * XYM(1) +
     &               XYN(2) * XYM(2) +
     &               XYN(3) * XYM(3)) / dsqrt(d1 * d2)

      Return
      End Function angle2Faces

C ================================================================
      Subroutine VecMul(a, b, c)
C ================================================================
      Real*8 a(3), b(3), c(3)

      c(1) = a(2) * b(3) - b(2) * a(3)
      c(2) = b(1) * a(3) - a(1) * b(3)
      c(3) = a(1) * b(2) - b(1) * a(2)
      Return
      End Subroutine VecMul

C ================================================================
      Real*8 Function DotMul(a, b)
C ================================================================
      Real*8 a(3), b(3)

      DotMul = a(1) * b(1) + a(2) * b(2) + a(3) * b(3)
      Return
      End Function DotMul
C
      End Module mba3d_auxSF


! { dg-do run }

program main
  implicit none
  real a, b
  real c
  !$acc declare create (c)

  a = 2.0
  b = 0.0

  !$acc parallel copy (a) create (b) default (none)
    b = a
    a = 1.0
    a = a + b
  !$acc end parallel

  if (a .ne. 3.0) call abort

  !$acc kernels copy (a) create (b) default (none)
    b = a
    a = 1.0
    a = a + b
  !$acc end kernels

  if (a .ne. 4.0) call abort

  !$acc parallel default (none) copy (a) create (b)
    b = a
    a = 1.0
    a = a + b
  !$acc end parallel

  if (a .ne. 5.0) call abort

  !$acc parallel default (none) copy (a)
    c = a
    a = 1.0
    a = a + c
  !$acc end parallel

  if (a .ne. 6.0) call abort

  !$acc data copy (a)
  !$acc parallel default (none)
    c = a
    a = 1.0
    a = a + c
  !$acc end parallel
  !$acc end data

  if (a .ne. 7.0) call abort

end program main


text:
!
! Empty Module for gcc compiler option
!

MODULE IFPORT
implicit none
END MODULE

max_stars_repo_name: solepomies/MAOOAM

max_issues_repo_name: solepomies/MAOOAM

max_forks_repo_name: solepomies/MAOOAM


text:
program hello
  implicit none
  write(*,*) 'Hello world!'
end program hello

max_stars_repo_name: xOluwaseyi/Hello-world

max_issues_repo_name: xOluwaseyi/Hello-world

max_forks_repo_name: xOluwaseyi/Hello-world


text:
      SUBROUTINE scle_velo()

      implicit none

      include 'p_input.inc'
      include 'p_array.inc'
      include 'p_const.inc'
      include 'cm_atom.inc'
      include 'cm_temp.inc'

      integer i,iatm

      real*8 instanke,xscale

      instanke=0.d0

      do i=1,natms

       iatm=atmtype(i)
       instanke=instanke
     $          +typmass(iatm)*((vvv(i,1)**2+vvv(i,2)**2+vvv(i,3)**2))

      enddo

      instanke=0.5*instanke
      xscale=sqrt(targetke/instanke)

      do i=1,natms

       vvv(i,1)=xscale*vvv(i,1)
       vvv(i,2)=xscale*vvv(i,2)
       vvv(i,3)=xscale*vvv(i,3)

      enddo

      return

      END
c $Id$

max_stars_repo_name: dinisAbranches/nwchem

max_issues_repo_name: dinisAbranches/nwchem

max_forks_repo_name: dinisAbranches/nwchem


text:
great(d1,k1).
great(ff1,ll1).
great(r1,bb1).
great(ii1,m1).
great(u1,p1).
great(hh1,l1).
great(g1,o1).
great(ii1,o1).
great(z1,j1).
great(y1,w1).
great(ff1,bb1).
great(cc1,s1).
great(j1,cc1).
great(hh1,ee1).
great(cc1,ll1).
great(u1,x1).
great(u1,v1).
great(f1,ee1).
great(d1,ff1).
great(hh1,k1).
great(d1,t1).
great(aa1,n1).
great(jj1,x1).
great(k1,cc1).
great(ee1,aa1).
great(z1,g1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!> \brief This module provides a subroutine to assemble heat capacity matrix
module m_heat_mat_ass_capacity
contains

  subroutine heat_mat_ass_capacity(hecMESH, hecMAT, fstrHEAT, delta_time)
    use m_fstr
    use m_heat_LIB_CAPACITY
    implicit none
    type(fstr_heat)          :: fstrHEAT
    type(hecmwST_matrix)     :: hecMAT
    type(hecmwST_local_mesh) :: hecMESH
    integer(kind=kint) :: i, in, j, nodLOCAL(20), ip, inod
    integer(kind=kint) :: itype, iS, iE, ic_type, icel, isect, IMAT, in0, nn
    real(kind=kreal) :: temp(20), lumped(120), mass(20*6, 20*6), ecoord(3,20)
    real(kind=kreal) :: delta_time, surf, THICK, ALFA, BETA

    beta = fstrHEAT%beta

    do itype = 1, hecMESH%n_elem_type
      iS = hecMESH%elem_type_index(itype-1) + 1
      iE = hecMESH%elem_type_index(itype  )
      ic_type = hecMESH%elem_type_item(itype)

      if (hecmw_is_etype_link(ic_type)) cycle
      if (hecmw_is_etype_patch(ic_type)) cycle
      if(ic_type == 3414) cycle

      !$omp parallel default(none), &
        !$omp&  private(icel,isect,IMAT,nn,temp,in0,i,j,nodLOCAL,ecoord,in,thick,surf,inod,lumped,mass), &
        !$omp&  shared(iS,iE,hecMESH,ic_type,hecMAT,fstrHEAT,delta_time)
      !$omp do
      do icel = iS, iE
        isect = hecMESH%section_ID(icel)
        IMAT = hecMESH%section%sect_mat_ID_item(isect)
        in0 = hecMESH%elem_node_index(icel-1)
        nn = hecmw_get_max_node(ic_type)
        do i = 1, nn
          nodLOCAL(i) = hecMESH%elem_node_item(in0+i)
          do j = 1, 3
            ecoord(j,i) = hecMESH%node(3*(nodLOCAL(i)-1)+j)
          enddo
          temp(i) = fstrHEAT%TEMP0(nodLOCAL(i))
        enddo

        lumped = 0.0d0
        if(ic_type == 231 .or. ic_type == 232 .or. ic_type == 241 .or. ic_type == 242)then
          in = hecMesh%section%sect_R_index(isect)
          thick = hecMESH%section%sect_R_item(in)
          call heat_capacity_C2(ic_type, nn, ecoord(1:2,1:nn), temp, IMAT, thick, lumped, mass, &
            fstrHEAT%CPtab(IMAT), fstrHEAT%CPtemp(IMAT,:), fstrHEAT%CPfuncA(IMAT,:) ,fstrHEAT%CPfuncB(IMAT,:), &
            fstrHEAT%RHOtab(IMAT), fstrHEAT%RHOtemp(IMAT,:), fstrHEAT%RHOfuncA(IMAT,:), fstrHEAT%RHOfuncB(IMAT,:))

        elseif(ic_type == 341 .or. ic_type == 342 .or. ic_type == 351 .or. ic_type == 352 .or. &
             & ic_type == 361 .or. ic_type == 362)then
          call heat_capacity_C3(ic_type, nn, ecoord(1:3,1:nn), temp, IMAT, lumped, mass, &
            fstrHEAT%CPtab(IMAT), fstrHEAT%CPtemp(IMAT,:), fstrHEAT%CPfuncA(IMAT,:) ,fstrHEAT%CPfuncB(IMAT,:), &
            fstrHEAT%RHOtab(IMAT), fstrHEAT%RHOtemp(IMAT,:), fstrHEAT%RHOfuncA(IMAT,:), fstrHEAT%RHOfuncB(IMAT,:))

        elseif (ic_type == 541) then
          ! skip warning

        elseif(ic_type == 731)then
          in = hecMesh%section%sect_R_index(isect)
          thick = hecMESH%section%sect_R_item(in)
          call heat_capacity_shell_731(ic_type, nn, ecoord(1:3,1:nn), temp, IMAT, thick, lumped, mass, &
            fstrHEAT%CPtab(IMAT), fstrHEAT%CPtemp(IMAT,:), fstrHEAT%CPfuncA(IMAT,:) ,fstrHEAT%CPfuncB(IMAT,:), &
            fstrHEAT%RHOtab(IMAT), fstrHEAT%RHOtemp(IMAT,:), fstrHEAT%RHOfuncA(IMAT,:), fstrHEAT%RHOfuncB(IMAT,:))

        elseif(ic_type == 741)then
          in = hecMesh%section%sect_R_index(isect)
          thick = hecMESH%section%sect_R_item(in)
          call heat_capacity_shell_741(ic_type, nn, ecoord(1:3,1:nn), temp, IMAT, thick, lumped, mass, &
            fstrHEAT%CPtab(IMAT), fstrHEAT%CPtemp(IMAT,:), fstrHEAT%CPfuncA(IMAT,:) ,fstrHEAT%CPfuncB(IMAT,:), &
            fstrHEAT%RHOtab(IMAT), fstrHEAT%RHOtemp(IMAT,:), fstrHEAT%RHOfuncA(IMAT,:), fstrHEAT%RHOfuncB(IMAT,:))

        elseif(ic_type == 111 .or. ic_type == 301 .or. ic_type == 611 .or. ic_type == 641)then
          in = hecMesh%section%sect_R_index(isect)
          surf = hecMESH%section%sect_R_item(in)
          call heat_capacity_C1(ic_type, nn, ecoord(1:3,1:nn), temp, IMAT, surf, lumped, mass, &
            fstrHEAT%CPtab(IMAT), fstrHEAT%CPtemp(IMAT,:), fstrHEAT%CPfuncA(IMAT,:) ,fstrHEAT%CPfuncB(IMAT,:), &
            fstrHEAT%RHOtab(IMAT), fstrHEAT%RHOtemp(IMAT,:), fstrHEAT%RHOfuncA(IMAT,:), fstrHEAT%RHOfuncB(IMAT,:))

        else
          write(*,*)"** error setMASS"
        endif

        do ip = 1, nn
          inod = nodLOCAL(ip)
          !$omp atomic
          hecMAT%D(inod) = hecMAT%D(inod) + lumped(ip) / delta_time
          !$omp atomic
          hecMAT%B(inod) = hecMAT%B(inod) + lumped(ip)*temp(ip) / delta_time
        enddo
      enddo
      !$omp end do
      !$omp end parallel
    enddo
  end subroutine heat_mat_ass_capacity
end module m_heat_mat_ass_capacity

max_stars_repo_name: fujitsu/FrontISTR

max_issues_repo_name: fujitsu/FrontISTR

max_forks_repo_name: fujitsu/FrontISTR


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                                !!
!!  This file is part of SciFT project                                               !!
!!  Copyright (c) 2011-2013 Nestor F. Aguirre (nfaguirrec@gmail.com)                 !!
!!                                                                                !!
!!  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 holders nor the names of its contributors   !!
!!     may be used to endorse or promote products derived from this software         !!
!!     without specific prior written permission.                                    !!
!!                                                                                !!
!!  The copyright holders provide no reassurances that the source code provided      !!
!!  does not infringe any patent, copyright, or any other intellectual property      !!
!!  rights of third parties.  The copyright holders disclaim any liability to any    !!
!!  recipient for claims brought against recipient by any third party for            !!
!!  infringement of that parties intellectual property rights.                       !!
!!                                                                                !!
!!  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 OWNER 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.                     !!
!!                                                                                !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! module Unit_
! 	implicit none
! 	private
!
! 	type, public :: Unit
! 		character(*) :: name
! 		character(*) :: symbol
! 		character(*) :: category
! 		real(8) :: value
! 	end type Unit
!
! end module Unit_

! Esto debe ser una clase singleton
! y entre lo que contendrá será
! el sistema de unidades con el que
! actualmente se está trabajando
! module SystemOfUnits_
! 	implicit none
! 	private
!
! 	type, public :: SystemOfUnits
! 		real(8) :: mass
! 	end type SystemOfUnits
!
! end module SystemOfUnits_

! module CM1Units_
! 	implicit none
! 	private
!
! 	!!
! 	! Mass   => amu
! 	! Lenght => Å
! 	! Time   => 9.14294598 10⁻¹¹ s
! 	! Energy => cm⁻¹
! 	!!
! 	real(8), public, parameter :: amu = 1.0_8
! 	real(8), public, parameter :: angs = 1.0_8
! 	real(8), public, parameter :: t_cm1 = 1.0_8
! 	real(8), public, parameter :: cm1 = 1.0_8
! 	real(8), public, parameter :: g = 1.0_8/1.660538782e-24
! 	real(8), public, parameter :: kg = 1.0_8/1.660538782e-27
! 	real(8), public, parameter :: m_au = 1.0_8/1822.88853_8
! 	real(8), public, parameter :: e_au = 1.0_8/4.55633538e-6
! 	real(8), public, parameter :: Eh = e_au
! 	real(8), public, parameter :: Ha = e_au
! end module CM1Units_

!!! Será renombrada UnitsManager
module UnitsConverter_
	use GOptions_
	use String_
! 	use StringRealPair_
! 	use StringRealMap_
	implicit none
	private

	public :: &
		sUnit, &
		UnitsConverter_test

	! -----------------
	! Esto no esta implementado
	! -----------------
	! Mass   => amu
	! Lenght => Å
	! Time   => 9.14294598 10⁻¹¹ s
	! Energy => cm⁻¹

	! -----------------
	! Atomic units
	! -----------------
	! Mass   => m_au
	! Lenght => bohr
	! Time   => t_au
	! Energy => Ha

! 	real(8), private :: LENGHT = 1.0_8
! 	real(8), private :: MASS = 1.0_8
! 	real(8), private :: ACTION = 1.0_8
! 	real(8), private :: CHARGE = 1.0_8

	integer, private, parameter :: UnitsConverter_NUNITS = 19

	character(10), private, parameter :: UnitsConverter_NAME(UnitsConverter_NUNITS) = [ &
		"amu", &
		"angs", &
		"eV", &
		"cm1", &
		"ryd", &
		"kcalMol", &
		"kJMol", &
		"bohr", &
		"g", &
		"kg", &
		"au", &
		"Eh", &
		"Ha", &
		"deg", &
		"kelvin", &
		"fs", &
		"ps", &
		"nm", &
		"microm" &
	]

	real(8), private, parameter :: UnitsConverter_VALUE(UnitsConverter_NUNITS) = [ &
		1822.88853_8, &                            ! amu
		1.0_8/0.52917726_8, &                      ! angs
		1.0_8/27.211396132_8, &                    ! eV
		1.0_8/219474.63068_8, &                    ! cm1
		1.0_8/2.0_8, &                             ! ryd
		1.0_8/627.50956_8, &                       ! kcalMol
		1.0_8/2625.49962_8, &                      ! kJMol
		1.0_8, &                                   ! bohr
		1.0_8/1.660538782e-24, &                   ! g    <<<<<<<<<<< Verificar
		1.0_8/1.660538782e-27, &                   ! kg   <<<<<<<<<<< Verificar
		1.0_8, &                                   ! au
		1.0_8, &                                   ! Eh
		1.0_8, &                                   ! Ha
		acos(-1.0_8)/180.0_8, &                    ! deg
		1.0/3.1577464d5, &                         ! kelvin
! 		1.0/0.02418884326505_8 &                   ! fs
		1.0d-15*1.0/2.418884326505d-17, &          ! fs
		1.0d-12*1.0/2.418884326505d-17, &          ! ps
		10.0_8/0.52917726_8, &                     ! nm
		1.0d3*10.0_8/0.52917726_8 &                ! microm
	]

	real(8), public, parameter :: amu       =   UnitsConverter_VALUE(1)
	real(8), public, parameter :: angs      =   UnitsConverter_VALUE(2)
	real(8), public, parameter :: eV        =   UnitsConverter_VALUE(3)
	real(8), public, parameter :: cm1       =   UnitsConverter_VALUE(4)
	real(8), public, parameter :: ryd       =   UnitsConverter_VALUE(5)
	real(8), public, parameter :: kcalMol   =   UnitsConverter_VALUE(6)
	real(8), public, parameter :: kJMol     =   UnitsConverter_VALUE(7)
	real(8), public, parameter :: bohr      =   UnitsConverter_VALUE(8)
	real(8), public, parameter :: g         =   UnitsConverter_VALUE(9)
	real(8), public, parameter :: kg        =   UnitsConverter_VALUE(10)
	real(8), public, parameter :: au        =   UnitsConverter_VALUE(11)
	real(8), public, parameter :: Eh        =   UnitsConverter_VALUE(12)
	real(8), public, parameter :: Ha        =   UnitsConverter_VALUE(13)
	real(8), public, parameter :: deg       =   UnitsConverter_VALUE(14)
	real(8), public, parameter :: kelvin    =   UnitsConverter_VALUE(15)
	real(8), public, parameter :: fs        =   UnitsConverter_VALUE(16)
	real(8), public, parameter :: ps        =   UnitsConverter_VALUE(17)
	real(8), public, parameter :: nm        =   UnitsConverter_VALUE(18)
	real(8), public, parameter :: microm    =   UnitsConverter_VALUE(19)

	contains

	!>
	!! @brief Constructor
	!!
	function sUnit( name ) result( output )
		character(*) :: name
		real(8) :: output

		integer :: i

		do i=1,UnitsConverter_NUNITS
			if ( trim(name) == trim(UnitsConverter_NAME(i)) ) then
				output = UnitsConverter_VALUE(i)
				return
			end if
		end do

		if( len_trim(name) > 0 ) call GOptions_warning( "Unknown unit ("//trim(name)//")", "UnitsConverter.sUnit()" )

		output = 1.0_8
	end function sUnit

	!>
	!! @brief Test method
	!!
	subroutine UnitsConverter_test()
! 		write(*,*) "         Lenght = ", 0.528e-8*cm
! 		write(*,*) "       Velocity = ", 2.18e8*cm/sec
		write(*,*) "         Energy = ", 27.21*eV
! 		write(*,*) "           Time = ", 24.2*asec
		write(*,*) "           Mass = ", 9.10938188e-31*kg
! 		write(*,*) "      Frequency = ", 4.13e16/sec
! 		write(*,*) " Electric field = ", 5.14e9*Volt/cm
! 		write(*,*) " Magnetic field = ", 2.35e5*Tesla
! 		write(*,*) "      Frequency = ", 6.6e-7*eV/GHz
! 		write(*,*) "Laser Intensity = ", 3.51e16*Wcm2
		write(*,*) ""
		write(*,*) "Lennard Jones units (argon)"
! 		write(*,*) "     Temperature = ", 120*K/LJ_T
! 		write(*,*) "Angular momentum = ", 0.02944*LJ_j
	end subroutine UnitsConverter_test

end module UnitsConverter_

max_stars_repo_name: nfaguirrec/scift


text:
      COMPLEX*16 FUNCTION C0_Wdec(p12,p22,Q2,m12,m22,m32)
*     ----------------------------------------------
* general C0; CalpHEP group, version 17/01/01
* Attention!  It doesn't work under physical treshold! To be improved!!!
*
      IMPLICIT NONE
      INTEGER I
      REAL*8 eps,p12,p22,Q2,p1dp2,DET3,SQD3,alpha,N,d,e,r0,pi,zet2
      REAL*8 rm12,rm22,rm32,rmb2,YR,RLOG,TMI,YPL,YMI,YDE,SIG
      REAL*8 DABS,DSQRT,DLOG,DDILOG_Wdec,DREAL,DIMAG
      COMPLEX*16 m12,m22,m32,mb2,ratio,LAMBDA,ALS,BETA
      COMPLEX*16 SQS,SQR,X0,X1,X2,X3,Y0,Y1,Y2,Y3,Y4,Y5,Y6,XSPENZ_Wdec
      COMPLEX*16 R21,R31,R13,R23,BETA21,BETA31,BETA13,BETA23,BETAZ
      COMPLEX*16 IMI,RKERN_Wdec,SKERN,SKERN1,SKERN2
      COMPLEX*16 AP,BP,CP,DP,SQDP,YP,YM,EP,DE,TERM,VELTET_Wdec
      COMPLEX*16 SQRT,LOG,DCMPLX
*
      DIMENSION AP(0:2),BP(0:2),CP(0:2),DP(0:2),SQDP(0:2)
      DIMENSION EP(0:2),DE(0:2),YR(0:2),YP(0:2),YM(0:2)
      DIMENSION SKERN1(0:2),SKERN2(0:2),SKERN(0:2)
*
      PARAMETER (eps=1D-20)
      PARAMETER (pi=3.14159265358979D0)
      PARAMETER (zet2=1.64493406684822618D0)
*
      rm12=DREAL(m12)
      rm22=DREAL(m22)
      rm32=DREAL(m32)
*
* Special cases, moved from C0_sp
*
*     c0(0,-rem2,Qs;wm2,em2,0)
*     c0(-rem2,0,Qs;0,wm2,em2)
*     c0(-rum2,0,Qs;0,um2,wm2) ???
*
*     print *,'args=',p12,p22,Q2,m12,m22,m32
      IF((p12.EQ.0D0.AND.-p22.EQ.rm22.AND.rm32.EQ.0D0).OR.
     &   (-p12.EQ.rm22.AND.p22.EQ.0D0.AND.rm12.EQ.0D0)) THEN
       IF(p12.EQ.0D0.AND.rm32.EQ.0D0) THEN
        mb2=m12
       ELSEIF(p22.EQ.0D0.AND.rm12.EQ.0D0) THEN
        mb2=m32
       ENDIF
       IF(rm22.LT.1D-2) THEN
        ratio=mb2/(Q2+mb2)
        C0_Wdec=1D0/Q2*(-1D0/2*LOG(ratio)**2
     &             -LOG(Q2/m22)*LOG(ratio)-XSPENZ_Wdec(ratio)+zet2)
        RETURN
       ELSE
        C0_Wdec=1D0/(Q2+m22)*(XSPENZ_Wdec(-Q2/m22*(mb2-m22)/(Q2+mb2))
     &    -XSPENZ_Wdec(-Q2/m22)-XSPENZ_Wdec((mb2-m22)/(Q2+mb2))+zet2)
        RETURN
       ENDIF
*
*     c0(-rem2,-rem2,Qs;cz2,em2,cz2)
*     c0(-rem2,-rem2,Qs;zm2,em2,cz2)
*     c0(-rem2,-rem2,Qs;cz2,em2,zm2)
*
      ELSEIF(p12.EQ.p22.AND.-p12.EQ.rm22.AND.(rm12.EQ.0D0.OR.rm32.EQ.0D0
     &                                       )) THEN
       IF(rm12.NE.0D0) mb2=m12
       IF(rm32.NE.0D0) mb2=m32
       IF(rm12.EQ.0D0.AND.rm32.EQ.0D0) mb2=m32
       rmb2=DREAL(mb2)
       IF(rmb2.LT.1D-2.AND.rm22.LT.1D-2) THEN
* GG light mass case ( LOG(Q2/M2)=LOG(S/M2)+SIG*DCMPLX(0D0,1D0)*PI )
        SIG=DABS(Q2)/Q2
        ALS=LOG(-Q2/m22)
        C0_Wdec=1D0/Q2*(1D0/2*ALS**2+SIG*DCMPLX(0D0,1D0)*PI*ALS+zet2)
        RETURN
       ELSEIF(rmb2.LT.1D-2.AND.rm22.GE.1D-2) THEN
* GG massive case
        BETA=SQRT(1D0+4D0*m22/Q2)
        Y1=-Q2/(2D0*m22)*(1D0+BETA)
        Y2=-Q2/(2D0*m22)*(1D0-BETA)
        C0_Wdec=-1D0/Q2/BETA*(XSPENZ_Wdec(1D0/Y1)-XSPENZ_Wdec(1D0/Y2)
     &                   +1D0/2*(LOG(-Y2)**2-LOG(-Y1)**2)
     &                   +2D0*(LOG(-Y1)-LOG(Y1))*LOG(Y2))
        RETURN
       ELSEIF(rmb2.GE.1D-2.AND.rm22.LT.1D-2) THEN
* ZG light mass case
        ALS=LOG(-Q2/m22)
        C0_Wdec=1D0/Q2*(ALS*LOG((Q2+mb2)/mb2)-zet2
     &         +XSPENZ_Wdec((mb2+Q2)/mb2))
        RETURN
       ELSEIF(rmb2.GE.1D-2.AND.rm22.GE.1D-2) THEN
* ZG massive case is computed below
        GO TO 100
       ELSE
        PRINT *,'Something wrong in c0_sp'
        STOP
       ENDIF
      ENDIF
*
 100  CONTINUE
*
* Cases p12=p22=0
*
      IF(p12.EQ.0D0.AND.p22.EQ.0D0.AND.rm12.EQ.rm32) THEN
*
* Case p12=p22=0 and m12,m22,m12
*
       SQR=SQRT(1D0+4D0*m12/Q2)
       X1=(1D0-SQR)/2D0
       X2=(1D0+SQR)/2D0
       X0=(m22-m12)/Q2
       X3=m22/(m22-m12)
       Y1=X1/X0
       Y2=(1D0-X1)/(1D0-X0)
       Y3=X2/X0
       Y4=(1D0-X2)/(1D0-X0)
       Y5=X3/X0
       Y6=(1D0-X3)/(1D0-X0)
*
       C0_Wdec=1D0/Q2*(XSPENZ_Wdec(1D0/(1D0-Y1))
     &           -XSPENZ_Wdec(1D0/(1D0-Y2))
     &           +XSPENZ_Wdec(1D0/(1D0-Y3))-XSPENZ_Wdec(1D0/(1D0-Y4))
     &           -XSPENZ_Wdec(1D0/(1D0-Y5))+XSPENZ_Wdec(1D0/(1D0-Y6)))
       RETURN
*
      ELSEIF(p12.EQ.0D0.AND.p22.EQ.0D0) THEN
*
* Case p12=p22=0 and arbitrary m12,m22,m32
*
       LAMBDA=Q2**2+2D0*Q2*(m12+m32)+(m12-m32)**2
       SQR=SQRT(LAMBDA)
       SQS=Q2+m12-m32
       X0=1D0+(m12-m22)/Q2
       X1=(SQS-SQR)/(2D0*Q2)
       X2=(SQS+SQR)/(2D0*Q2)
       X3=m32/(m32-m22)
       C0_Wdec=1D0/Q2*(XSPENZ_Wdec((X0-1D0)/(X0-X1))
     &           -XSPENZ_Wdec(X0/(X0-X1))
     &           +XSPENZ_Wdec((X0-1D0)/(X0-X2))-XSPENZ_Wdec(X0/(X0-X2))
     &           -XSPENZ_Wdec((X0-1D0)/(X0-X3))+XSPENZ_Wdec(X0/(X0-X3)))
       RETURN
      ENDIF
*
* Pinch of WW-direct box, cat0wb = C0WWB(0d0,-rtm2,-t,cz2,wm2,bm2)
*                         cat0wb = C0WWB(-rtm2,0d0,-t,bm2,wm2,cz2)
*                         cat0w0 = C0WWB(0d0,-rtm2,-t,cz2,wm2,cz2)
*
       IF((p12.EQ.0D0.AND.-p22.NE.rm32.AND.Q2.NE.0D0.AND.rm12.EQ.0D0)
     &.OR.(p22.EQ.0D0.AND.-p12.NE.rm12.AND.Q2.NE.0D0.AND.rm32.EQ.0D0)
     &   ) THEN
       IF(p12.EQ.0D0.AND.rm12.EQ.0D0) THEN
        TMI=p22-Q2
        YPL=(rm22-p22-rm32+DSQRT((rm22-p22-rm32)**2+4D0*rm22*p22))
     &                                                 /2D0/(-p22)
        YMI=rm22/(-p22)/YPL
        YDE=1D0+rm32/Q2
       ENDIF
       IF(p22.EQ.0D0.AND.rm32.EQ.0D0) THEN
        TMI=p12-Q2
        YPL=(rm22-p12-rm12+DSQRT((rm22-p12-rm12)**2+4D0*rm22*p12))
     &                                                 /2D0/(-p12)
        YMI=rm22/(-p12)/YPL
        YDE=1D0+rm12/Q2
       ENDIF
       C0_Wdec=DCMPLX(1D0/TMI*(-DDILOG_Wdec((rm22+TMI)/(rm22+TMI*YPL))
     &                    +DDILOG_Wdec((rm22    )/(rm22+TMI*YPL))
     &                    -DDILOG_Wdec((rm22+TMI)/(rm22+TMI*YMI))
     &                    +DDILOG_Wdec((rm22    )/(rm22+TMI*YMI))
     &                    +DDILOG_Wdec((rm22+TMI)/(rm22+TMI*YDE))
     &                    -DDILOG_Wdec((rm22    )/(rm22+TMI*YDE))
     &                    +DDILOG_Wdec(1D0+TMI/rm22)-zet2),0D0  )
       RETURN
*
* Two W-decay cases
*
      ELSEIF(-p12.EQ.rm22.AND.p22.EQ.0D0.AND.-Q2.EQ.rm32) THEN
*
       R21=m22/m12
       R31=m32/m12
       BETA21=SQRT(1D0-4D0*R21)
       BETA31=SQRT(1D0-4D0*R31)
       C0_Wdec=1D0/(m32-m22)*(LOG((1D0+BETA31)/2D0)*LOG(2D0*R31
     &         /(1D0+BETA31))
     &                 -LOG((1D0+BETA21)/2D0)*LOG(2D0*R21/(1D0+BETA21)))
       RETURN
*
      ELSEIF(p12.EQ.0D0.AND.-p22.EQ.rm22.AND.-Q2.EQ.rm12) THEN
*
       R13=m12/m32
       R23=m22/m32
       BETA13=SQRT(1D0-4D0*R13)
       BETA23=SQRT(1D0-4D0*R23)
       C0_Wdec=1D0/(m12-m22)*(LOG((1D0+BETA13)/2D0)*LOG(2D0*R13
     &        /(1D0+BETA13))
     &                 -LOG((1D0+BETA23)/2D0)*LOG(2D0*R23/(1D0+BETA23)))
       RETURN
      ENDIF
*
* Variations of exact case
*
* 1) c0(0d0,-rtm2,-t,cz2,zm2,tm2)
*
       IF((p12.EQ.0D0.AND.-p22.EQ.rm32.AND.Q2.NE.0D0.AND.rm12.EQ.0D0)
     &.OR.(p12.EQ.0D0.AND.-p22.EQ.rm12.AND.Q2.NE.0D0.AND.rm32.EQ.0D0)
     &   ) THEN
*
* Case p12.eq.0.and.(-p22.eq.rm12.or.rm32).and.Q2.ne.0
*
       BETAZ=SQRT(1D0+4D0*p22/m22)
       Y0=m22/(Q2-p22)
       Y1=2D0/(1D0-BETAZ)
       Y2=2D0/(1D0+BETAZ)
       C0_Wdec=1D0/(Q2-p22)*(zet2-XSPENZ_Wdec(1D0-1D0/Y0)
     &                     -XSPENZ_Wdec((1D0-Y0)*Q2/(-p22+Q2-Y0*Q2))
     &                     +XSPENZ_Wdec((   -Y0)*Q2/(-p22+Q2-Y0*Q2))
     &  +XSPENZ_Wdec((1D0-Y0)/(Y1-Y0))-XSPENZ_Wdec((   -Y0)/(Y1-Y0))
     &  +XSPENZ_Wdec((1D0-Y0)/(Y2-Y0))-XSPENZ_Wdec((   -Y0)/(Y2-Y0)))
       RETURN
*
* 2) c0(-rtm2,0d0,-t,tm2,zm2,cz2)
*
      ELSEIF((-p12.EQ.rm12.AND.p22.EQ.0D0.AND.Q2.NE.0D0.AND.rm32.EQ.0D0)
     &   .OR.(-p12.EQ.rm32.AND.p22.EQ.0D0.AND.Q2.NE.0D0.AND.rm12.EQ.0D0)
     &      ) THEN
*
* Case (-p12.eq.rm12.or.rm32).and.p22.eq.0.and.Q2.ne.0
*
       BETAZ=SQRT(1D0+4D0*p12/m22)
       Y0=m22/(Q2-p12)
       Y1=2D0/(1D0-BETAZ)
       Y2=2D0/(1D0+BETAZ)
       C0_Wdec=1D0/(Q2-p12)*(zet2-XSPENZ_Wdec(1D0-1D0/Y0)
     &                     -XSPENZ_Wdec((1D0-Y0)*Q2/(-p12+Q2-Y0*Q2))
     &                     +XSPENZ_Wdec((   -Y0)*Q2/(-p12+Q2-Y0*Q2))
     &  +XSPENZ_Wdec((1D0-Y0)/(Y1-Y0))-XSPENZ_Wdec((   -Y0)/(Y1-Y0))
     &  +XSPENZ_Wdec((1D0-Y0)/(Y2-Y0))-XSPENZ_Wdec((   -Y0)/(Y2-Y0)))
       RETURN
* 3)
      ELSEIF(p12.NE.0D0.AND.p22.EQ.0D0.AND.-Q2.EQ.rm12.AND.rm22.EQ.0D0
     &      ) THEN
*
* Case p12.ne.0.and.-Q2.eq.rm12.and.p2.eq.0
*
       BETAZ=SQRT(1D0+4D0*Q2/m32)
       Y0=m32/(p12-Q2)
       Y1=2D0/(1D0-BETAZ)
       Y2=2D0/(1D0+BETAZ)
       C0_Wdec=1D0/(p12-Q2)*(zet2-XSPENZ_Wdec(1D0-1D0/Y0)
     &                     -XSPENZ_Wdec((1D0-Y0)*p12/(-Q2+p12-Y0*p12))
     &                     +XSPENZ_Wdec((   -Y0)*p12/(-Q2+p12-Y0*p12))
     &  +XSPENZ_Wdec((1D0-Y0)/(Y1-Y0))-XSPENZ_Wdec((   -Y0)/(Y1-Y0))
     &  +XSPENZ_Wdec((1D0-Y0)/(Y2-Y0))-XSPENZ_Wdec((   -Y0)/(Y2-Y0)))
       RETURN
* 4)
      ELSEIF(p12.EQ.0D0.AND.p22.NE.0D0.AND.-Q2.EQ.rm32.AND.rm22.EQ.0D0
     &      ) THEN
*
* Case -p12.eq.0.and.p22.ne.0.and.Q2.eq.rm32; PROBLEMS HERE
*
       BETAZ=SQRT(1D0+4D0*Q2/m12)
       Y0=m12/(p22-Q2)
       Y1=2D0/(1D0-BETAZ)
       Y2=2D0/(1D0+BETAZ)
       C0_Wdec=1D0/(p22-Q2)*(zet2-XSPENZ_Wdec(1D0-1D0/Y0)
     &                     -XSPENZ_Wdec((1D0-Y0)*p22/(-Q2+p22-Y0*p22))
     &                     +XSPENZ_Wdec((   -Y0)*p22/(-Q2+p22-Y0*p22))
     &  +XSPENZ_Wdec((1D0-Y0)/(Y1-Y0))-XSPENZ_Wdec((   -Y0)/(Y1-Y0))
     &  +XSPENZ_Wdec((1D0-Y0)/(Y2-Y0))-XSPENZ_Wdec((   -Y0)/(Y2-Y0)))
       RETURN
* 5)
      ELSEIF(p12.NE.0D0.AND.p22.NE.0D0) THEN
*
* Case p12.ne.0.and.p22.ne.0; or GENERAL case
*
      IMI =DCMPLX(0D0,1D0)
      p1dp2=(Q2-p12-p22)/2D0
      DET3 =p1dp2**2-p12*p22
      SQD3 =DSQRT(DET3)
      IF(p1dp2.ge.0d0) THEN
       alpha=(p1dp2+SQD3)/(-p12)
      ELSE
       alpha=(-p22)/(p1dp2-SQD3)
      ENDIF
      N=2D0*SQD3
c     print *,'alpha,N=',alpha,N
      d=DREAL(p22+m22-m32)
      e=DREAL(Q2-p22+m12-m22)
      r0=(-d-e*alpha)/N
*
      YR(0) =-r0/alpha
      YR(1) = r0/(1d0-alpha)
      YR(2) = r0+alpha
c     print *,'Roots=',YR(0),YR(1),YR(2)
*
      AP(0) =-p22
      BP(0) = p22+m22-m32
      CP(0) = m32
*
      AP(1) =-Q2
      BP(1) = Q2+m12-m32
      CP(1) = m32
*
      AP(2) =-p12
      BP(2) = p12+m12-m22
      CP(2) = m22
*
      EP(0)=eps
      EP(1)=eps
      EP(2)=eps
*
      DE(0)=-DIMAG(AP(0)*YR(0)**2+BP(0)*YR(0)+CP(0)-IMI*eps)
      DE(1)=-DIMAG(AP(1)*YR(1)**2+BP(1)*YR(1)+CP(1)-IMI*eps)
      DE(2)=-DIMAG(AP(2)*YR(2)**2+BP(2)*YR(2)+CP(2)-IMI*eps)
*
      C0_Wdec=0D0
*
      DO I=0,2,1
      DP(I)=BP(I)**2-4D0*AP(I)*CP(I)
c     print *,'I DP(I)=',I,DP(I)
      SQDP(I)=SQRT(DP(I))
      YP(I)=(-BP(I)+SQDP(I))/2D0/AP(I)
c     print *,'I YP(I)=',I,YP(I)
      IF(DABS(DREAL(YP(I))).LT.1D-7) THEN
       YP(I)=-2D0*CP(I)/(BP(I)+SQDP(I))
c      print *,'I YP(I)=',I,YP(I)
      ENDIF
      YM(I)=(-BP(I)-SQDP(I))/2D0/AP(I)
c     print *,'I YM(I)=',I,YM(I)
      IF(DABS(DREAL(YM(I))).LT.1D-7) THEN
       YM(I)=-2D0*CP(I)/(BP(I)-SQDP(I))
c      print *,'I YM(I)=',I,YM(I)
      ENDIF
      IF(Q2.GT.0D0) THEN
        RLOG=DLOG(DABS(1D0-1D0/YR(I)))
      ELSE
        RLOG=DLOG(DABS(1D0-1D0/YR(I)))
      ENDIF
      SKERN1(I)=-(+VELTET_Wdec(-YP(I),-YM(I))
     &         -VELTET_Wdec(YR(I)-YP(I),YR(I)-YM(I))
     &           )*RLOG
      SKERN2(I)=-(-VELTET_Wdec(AP(I)-IMI*EP(I),1D0/(AP(I)-IMI*DE(I))))
     &         *RLOG
      SKERN(I)=SKERN1(I)+SKERN2(I)
*
* Analyser SKERN1
*
*     IF(DREAL(SKERN1(I)).NE.0D0.OR.DIMAG(SKERN1(I)).NE.0D0)
*    &              PRINT *,'I,Re(SKERN1(I)),Im(SKERN1(I))=',I,SKERN1(I)
*
* Analyser SKERN2
*
*     IF(DREAL(SKERN2(I)).NE.0D0.OR.DIMAG(SKERN2(I)).NE.0D0)
*    &              PRINT *,'I,Re(SKERN2(I)),Im(SKERN2(I))=',I,SKERN2(I)
*
      TERM=(-1)**I/N*(RKERN_Wdec(YR(I),YP(I))+RKERN_Wdec(YR(I),YM(I))
     &    +SKERN(I))
      C0_Wdec=C0_Wdec+TERM
      ENDDO
*
      ELSE
       PRINT *,'Combination of arguments is not yet foreseen!'
       STOP
      ENDIF
*
      RETURN
      END

      COMPLEX*16 FUNCTION RKERN_Wdec(Y,YL)
*
      IMPLICIT NONE
      REAL*8 Y,DREAL,DIMAG
      COMPLEX*16 YL,RKERNR,RKERNI,XSPENZ_Wdec,LOG,VELTET_Wdec
*
      RKERNR=XSPENZ_Wdec(Y/(Y-YL))-XSPENZ_Wdec((Y-1D0)/(Y-YL))
      RKERNI=VELTET_Wdec(   -YL,1D0/(Y-YL))*LOG(Y/(Y-YL))
     &      -VELTET_Wdec(1D0-YL,1D0/(Y-YL))*LOG((Y-1D0)/(Y-YL))
      RKERN_Wdec=RKERNR+RKERNI
*
* Analyser RKERN
*
      IF(DREAL(RKERNI).NE.0D0.OR.DIMAG(RKERNI).NE.0D0)
     &                                   PRINT *,'Re(RI),Im(RI)=',RKERNI
*
      RETURN
      END

      COMPLEX*16 FUNCTION VELTET_Wdec(A,B)
*
      IMPLICIT NONE
      REAL*8 PI,IMA,IMB,IMAB,TET_Wdec

      COMPLEX*16 A,B,IMI,DCMPLX
*
      DATA PI/3.14159274101D0/
*
      IMA =DIMAG(A)
      IMB =DIMAG(B)
      IMAB=DIMAG(A*B)
      IMI =DCMPLX(0D0,1D0)
      VELTET_Wdec=2D0*IMI*PI
     &           *(TET_Wdec(-IMA)*TET_Wdec(-IMB)*TET_Wdec( IMAB)
     &                  -TET_Wdec( IMA)*TET_Wdec( IMB)*TET_Wdec(-IMAB))
*
      RETURN
      END

      REAL*8 FUNCTION TET_Wdec(X)
*
      IMPLICIT NONE
      REAL*8 X
*
      IF(X.LT.0D0) THEN
        TET_Wdec=0D0
      ELSE
        TET_Wdec=1D0
      ENDIF
*
      RETURN
      END

max_stars_repo_name: lukaszgaza/winhacplusplus

max_issues_repo_name: lukaszgaza/winhacplusplus

max_forks_repo_name: lukaszgaza/winhacplusplus


text:
      DOUBLE PRECISION FUNCTION iau_EECT00 ( DATE1, DATE2 )
*+
*  - - - - - - - - - - -
*   i a u _ E E C T 0 0
*  - - - - - - - - - - -
*
*  Equation of the equinoxes complementary terms, consistent with
*  IAU 2000 resolutions.
*
*  This routine is part of the International Astronomical Union's
*  SOFA (Standards of Fundamental Astronomy) software collection.
*
*  Status:  canonical model.
*
*  Given:
*     DATE1,DATE2   d    TT as a 2-part Julian Date (Note 1)
*
*  Returned:
*     iau_EECT00    d    complementary terms (Note 2)
*
*  Notes:
*
*  1) The TT date DATE1+DATE2 is a Julian Date, apportioned in any
*     convenient way between the two arguments.  For example,
*     JD(TT)=2450123.7 could be expressed in any of these ways,
*     among others:
*
*            DATE1          DATE2
*
*         2450123.7D0        0D0        (JD method)
*          2451545D0      -1421.3D0     (J2000 method)
*         2400000.5D0     50123.2D0     (MJD method)
*         2450123.5D0       0.2D0       (date & time method)
*
*     The JD method is the most natural and convenient to use in
*     cases where the loss of several decimal digits of resolution
*     is acceptable.  The J2000 method is best matched to the way
*     the argument is handled internally and will deliver the
*     optimum resolution.  The MJD method and the date & time methods
*     are both good compromises between resolution and convenience.
*
*  2) The "complementary terms" are part of the equation of the
*     equinoxes (EE), classically the difference between apparent and
*     mean Sidereal Time:
*
*        GAST = GMST + EE
*
*     with:
*
*        EE = dpsi * cos(eps)
*
*     where dpsi is the nutation in longitude and eps is the obliquity
*     of date.  However, if the rotation of the Earth were constant in
*     an inertial frame the classical formulation would lead to apparent
*     irregularities in the UT1 timescale traceable to side-effects of
*     precession-nutation.  In order to eliminate these effects from
*     UT1, "complementary terms" were introduced in 1994 (IAU, 1994) and
*     took effect from 1997 (Capitaine and Gontier, 1993):
*
*        GAST = GMST + CT + EE
*
*     By convention, the complementary terms are included as part of the
*     equation of the equinoxes rather than as part of the mean Sidereal
*     Time.  This slightly compromises the "geometrical" interpretation
*     of mean sidereal time but is otherwise inconsequential.
*
*     The present routine computes CT in the above expression,
*     compatible with IAU 2000 resolutions (Capitaine et al., 2002, and
*     IERS Conventions 2003).
*
*  Called:
*     iau_FAL03    mean anomaly of the Moon
*     iau_FALP03   mean anomaly of the Sun
*     iau_FAF03    mean argument of the latitude of the Moon
*     iau_FAD03    mean elongation of the Moon from the Sun
*     iau_FAOM03   mean longitude of the Moon's ascending node
*     iau_FAVE03   mean longitude of Venus
*     iau_FAE03    mean longitude of Earth
*     iau_FAPA03   general accumulated precession in longitude
*
*  References:
*
*     Capitaine, N. & Gontier, A.-M., Astron. Astrophys., 275,
*     645-650 (1993)
*
*     Capitaine, N., Wallace, P.T. and McCarthy, D.D., "Expressions to
*     implement the IAU 2000 definition of UT1", Astronomy &
*     Astrophysics, 406, 1135-1149 (2003)
*
*     IAU Resolution C7, Recommendation 3 (1994)
*
*     McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
*     IERS Technical Note No. 32, BKG (2004)
*
*  This revision:  2013 December 2
*
*  SOFA release 2017-04-20
*
*  Copyright (C) 2017 IAU SOFA Board.  See notes at end.
*
*-----------------------------------------------------------------------

      IMPLICIT NONE

      DOUBLE PRECISION DATE1, DATE2

*  Arcseconds to radians
      DOUBLE PRECISION DAS2R
      PARAMETER ( DAS2R = 4.848136811095359935899141D-6 )

*  Reference epoch (J2000.0), JD
      DOUBLE PRECISION DJ00
      PARAMETER ( DJ00 = 2451545D0 )

*  Days per Julian century
      DOUBLE PRECISION DJC
      PARAMETER ( DJC = 36525D0 )

*  Time since J2000.0, in Julian centuries
      DOUBLE PRECISION T

*  Miscellaneous
      INTEGER I, J
      DOUBLE PRECISION A, S0, S1
      DOUBLE PRECISION iau_FAL03, iau_FALP03, iau_FAF03,
     :                 iau_FAD03, iau_FAOM03, iau_FAVE03, iau_FAE03,
     :                 iau_FAPA03

*  Fundamental arguments
      DOUBLE PRECISION FA(14)

*  -----------------------------------------
*  The series for the EE complementary terms
*  -----------------------------------------

*  Number of terms in the series
      INTEGER NE0, NE1
      PARAMETER ( NE0=33, NE1=1 )

*  Coefficients of l,l',F,D,Om,LVe,LE,pA
      INTEGER KE0 ( 8, NE0 ),
     :        KE1 ( 8, NE1 )

*  Sine and cosine coefficients
      DOUBLE PRECISION SE0 ( 2, NE0 ),
     :                 SE1 ( 2, NE1 )

*  Argument coefficients for t^0
      DATA ( ( KE0(I,J), I=1,8), J=1,10 ) /
     :  0,  0,  0,  0,  1,  0,  0,  0,
     :  0,  0,  0,  0,  2,  0,  0,  0,
     :  0,  0,  2, -2,  3,  0,  0,  0,
     :  0,  0,  2, -2,  1,  0,  0,  0,
     :  0,  0,  2, -2,  2,  0,  0,  0,
     :  0,  0,  2,  0,  3,  0,  0,  0,
     :  0,  0,  2,  0,  1,  0,  0,  0,
     :  0,  0,  0,  0,  3,  0,  0,  0,
     :  0,  1,  0,  0,  1,  0,  0,  0,
     :  0,  1,  0,  0, -1,  0,  0,  0 /
      DATA ( ( KE0(I,J), I=1,8), J=11,20 ) /
     :  1,  0,  0,  0, -1,  0,  0,  0,
     :  1,  0,  0,  0,  1,  0,  0,  0,
     :  0,  1,  2, -2,  3,  0,  0,  0,
     :  0,  1,  2, -2,  1,  0,  0,  0,
     :  0,  0,  4, -4,  4,  0,  0,  0,
     :  0,  0,  1, -1,  1, -8, 12,  0,
     :  0,  0,  2,  0,  0,  0,  0,  0,
     :  0,  0,  2,  0,  2,  0,  0,  0,
     :  1,  0,  2,  0,  3,  0,  0,  0,
     :  1,  0,  2,  0,  1,  0,  0,  0 /
      DATA ( ( KE0(I,J), I=1,8), J=21,30 ) /
     :  0,  0,  2, -2,  0,  0,  0,  0,
     :  0,  1, -2,  2, -3,  0,  0,  0,
     :  0,  1, -2,  2, -1,  0,  0,  0,
     :  0,  0,  0,  0,  0,  8,-13, -1,
     :  0,  0,  0,  2,  0,  0,  0,  0,
     :  2,  0, -2,  0, -1,  0,  0,  0,
     :  1,  0,  0, -2,  1,  0,  0,  0,
     :  0,  1,  2, -2,  2,  0,  0,  0,
     :  1,  0,  0, -2, -1,  0,  0,  0,
     :  0,  0,  4, -2,  4,  0,  0,  0 /
      DATA ( ( KE0(I,J), I=1,8), J=31,NE0 ) /
     :  0,  0,  2, -2,  4,  0,  0,  0,
     :  1,  0, -2,  0, -3,  0,  0,  0,
     :  1,  0, -2,  0, -1,  0,  0,  0 /

*  Argument coefficients for t^1
      DATA ( ( KE1(I,J), I=1,8), J=1,NE1 ) /
     :  0,  0,  0,  0,  1,  0,  0,  0 /

*  Sine and cosine coefficients for t^0
      DATA ( ( SE0(I,J), I=1,2), J =    1,   10 ) /
     :            +2640.96D-6,          -0.39D-6,
     :              +63.52D-6,          -0.02D-6,
     :              +11.75D-6,          +0.01D-6,
     :              +11.21D-6,          +0.01D-6,
     :               -4.55D-6,          +0.00D-6,
     :               +2.02D-6,          +0.00D-6,
     :               +1.98D-6,          +0.00D-6,
     :               -1.72D-6,          +0.00D-6,
     :               -1.41D-6,          -0.01D-6,
     :               -1.26D-6,          -0.01D-6 /
      DATA ( ( SE0(I,J), I=1,2), J =   11,   20 ) /
     :               -0.63D-6,          +0.00D-6,
     :               -0.63D-6,          +0.00D-6,
     :               +0.46D-6,          +0.00D-6,
     :               +0.45D-6,          +0.00D-6,
     :               +0.36D-6,          +0.00D-6,
     :               -0.24D-6,          -0.12D-6,
     :               +0.32D-6,          +0.00D-6,
     :               +0.28D-6,          +0.00D-6,
     :               +0.27D-6,          +0.00D-6,
     :               +0.26D-6,          +0.00D-6 /
      DATA ( ( SE0(I,J), I=1,2), J =   21,   30 ) /
     :               -0.21D-6,          +0.00D-6,
     :               +0.19D-6,          +0.00D-6,
     :               +0.18D-6,          +0.00D-6,
     :               -0.10D-6,          +0.05D-6,
     :               +0.15D-6,          +0.00D-6,
     :               -0.14D-6,          +0.00D-6,
     :               +0.14D-6,          +0.00D-6,
     :               -0.14D-6,          +0.00D-6,
     :               +0.14D-6,          +0.00D-6,
     :               +0.13D-6,          +0.00D-6 /
      DATA ( ( SE0(I,J), I=1,2), J =   31,  NE0 ) /
     :               -0.11D-6,          +0.00D-6,
     :               +0.11D-6,          +0.00D-6,
     :               +0.11D-6,          +0.00D-6 /

*  Sine and cosine coefficients for t^1
      DATA ( ( SE1(I,J), I=1,2), J =    1,  NE1 ) /
     :               -0.87D-6,          +0.00D-6 /

* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

*  Interval between fundamental epoch J2000.0 and current date (JC).
      T = ( ( DATE1-DJ00 ) + DATE2 ) / DJC

*  Fundamental Arguments (from IERS Conventions 2003)

*  Mean anomaly of the Moon.
      FA(1) = iau_FAL03 ( T )

*  Mean anomaly of the Sun.
      FA(2) = iau_FALP03 ( T )

*  Mean longitude of the Moon minus that of the ascending node.
      FA(3) = iau_FAF03 ( T )

*  Mean elongation of the Moon from the Sun.
      FA(4) = iau_FAD03 ( T )

*  Mean longitude of the ascending node of the Moon.
      FA(5) = iau_FAOM03 ( T )

*  Mean longitude of Venus.
      FA(6) = iau_FAVE03 ( T )

*  Mean longitude of Earth.
      FA(7) = iau_FAE03 ( T )

*  General precession in longitude.
      FA(8) = iau_FAPA03 ( T )

*  Evaluate the EE complementary terms.
      S0 = 0D0
      S1 = 0D0

      DO 2 I = NE0,1,-1
         A = 0D0
         DO 1 J=1,8
            A = A + DBLE(KE0(J,I))*FA(J)
 1       CONTINUE
         S0 = S0 + ( SE0(1,I)*SIN(A) + SE0(2,I)*COS(A) )
 2    CONTINUE
      DO 4 I = NE1,1,-1
         A = 0D0
         DO 3 J=1,8
            A = A + DBLE(KE1(J,I))*FA(J)
 3       CONTINUE
         S1 = S1 + ( SE1(1,I)*SIN(A) + SE1(2,I)*COS(A) )
 4    CONTINUE
      iau_EECT00 = ( S0 + S1 * T ) * DAS2R

*  Finished.

*+----------------------------------------------------------------------
*
*  Copyright (C) 2017
*  Standards Of Fundamental Astronomy Board
*  of the International Astronomical Union.
*
*  =====================
*  SOFA Software License
*  =====================
*
*  NOTICE TO USER:
*
*  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING SIX TERMS AND
*  CONDITIONS WHICH APPLY TO ITS USE.
*
*  1. The Software is owned by the IAU SOFA Board ("SOFA").
*
*  2. Permission is granted to anyone to use the SOFA software for any
*     purpose, including commercial applications, free of charge and
*     without payment of royalties, subject to the conditions and
*     restrictions listed below.
*
*  3. You (the user) may copy and distribute SOFA source code to others,
*     and use and adapt its code and algorithms in your own software,
*     on a world-wide, royalty-free basis.  That portion of your
*     distribution that does not consist of intact and unchanged copies
*     of SOFA source code files is a "derived work" that must comply
*     with the following requirements:
*
*     a) Your work shall be marked or carry a statement that it
*        (i) uses routines and computations derived by you from
*        software provided by SOFA under license to you; and
*        (ii) does not itself constitute software provided by and/or
*        endorsed by SOFA.
*
*     b) The source code of your derived work must contain descriptions
*        of how the derived work is based upon, contains and/or differs
*        from the original SOFA software.
*
*     c) The names of all routines in your derived work shall not
*        include the prefix "iau" or "sofa" or trivial modifications
*        thereof such as changes of case.
*
*     d) The origin of the SOFA components of your derived work must
*        not be misrepresented;  you must not claim that you wrote the
*        original software, nor file a patent application for SOFA
*        software or algorithms embedded in the SOFA software.
*
*     e) These requirements must be reproduced intact in any source
*        distribution and shall apply to anyone to whom you have
*        granted a further right to modify the source code of your
*        derived work.
*
*     Note that, as originally distributed, the SOFA software is
*     intended to be a definitive implementation of the IAU standards,
*     and consequently third-party modifications are discouraged.  All
*     variations, no matter how minor, must be explicitly marked as
*     such, as explained above.
*
*  4. You shall not cause the SOFA software to be brought into
*     disrepute, either by misuse, or use for inappropriate tasks, or
*     by inappropriate modification.
*
*  5. The SOFA software is provided "as is" and SOFA makes no warranty
*     as to its use or performance.   SOFA does not and cannot warrant
*     the performance or results which the user may obtain by using the
*     SOFA software.  SOFA makes no warranties, express or implied, as
*     to non-infringement of third party rights, merchantability, or
*     fitness for any particular purpose.  In no event will SOFA be
*     liable to the user for any consequential, incidental, or special
*     damages, including any lost profits or lost savings, even if a
*     SOFA representative has been advised of such damages, or for any
*     claim by any third party.
*
*  6. The provision of any version of the SOFA software under the terms
*     and conditions specified herein does not imply that future
*     versions will also be made available under the same terms and
*     conditions.
*
*  In any published work or commercial product which uses the SOFA
*  software directly, acknowledgement (see www.iausofa.org) is
*  appreciated.
*
*  Correspondence concerning SOFA software should be addressed as
*  follows:
*
*      By email:  sofa@ukho.gov.uk
*      By post:   IAU SOFA Center
*                 HM Nautical Almanac Office
*                 UK Hydrographic Office
*                 Admiralty Way, Taunton
*                 Somerset, TA1 2DN
*                 United Kingdom
*
*-----------------------------------------------------------------------

      END

ext: for


text:
**==cther.spg  processed by SPAG 4.50J  at 14:49 on 30 Jun 1995
      FUNCTION CTHER(N,J,T)
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL a , cr , CTHER , T , t4
      INTEGER in , J , N
C*** End of declarations inserted by SPAG
C  FITS TO SCOTT'S RATES FOR  HE I + J   TO   HE II + J-1
C  1000  TO 100,000 K THROUGH ARGON
      DIMENSION a(2,9)
      DATA a/4*0. , .059 , .00001 , 1.1 , 0.7 , .00001 , 1.7 , .075 ,
     &     2.2 , .096 , 1.2 , 1.1 , .0008 , .00001 , 1.0/
C
      CTHER = 0.
      IF ( J.EQ.1 ) RETURN
      t4 = AMAX1(.0001*T,0.1)
      t4 = AMIN1(t4,10.)
      cr = 5.4E-10*t4
      IF ( J.EQ.2 .OR. N.GT.18 ) RETURN
      CTHER = AMAX1(cr,(-.1+.4*J)*1.0E-9)
      IF ( J.GE.6 ) RETURN
      IF ( J.EQ.3 ) THEN
         CTHER = 0.
C  INCLUDE ALBERT'S RATE TO 1D OF N+
         IF ( N.EQ.7 ) CTHER = 3.0E-10*(1.+.26*t4) + 6.2E-11
         IF ( N.EQ.8 ) CTHER = 3.3E-10*t4**.7
         IF ( N.EQ.18 ) CTHER = 1.3E-10
         RETURN
      ELSEIF ( N.NE.7 ) THEN
         in = N/2
         CTHER = a(J-3,in)*1.0E-9
         IF ( J.NE.5 ) THEN
            IF ( N.EQ.6 ) CTHER = 5.1E-11*t4**1.46
            IF ( N.EQ.14 ) CTHER = 9.6E-10*t4**.55
            IF ( N.EQ.16 ) CTHER = 1.1E-9*t4**.63
            RETURN
         ENDIF
      ELSE
         CTHER = 1.1E-10
         IF ( J.EQ.5 ) CTHER = 2.0E-9
         RETURN
      ENDIF
      IF ( N.EQ.18 ) CTHER = 1.0E-9*t4**.91
      RETURN
      END

max_stars_repo_name: DougBurke/xspeclmodels

max_issues_repo_name: DougBurke/xspeclmodels

max_forks_repo_name: DougBurke/xspeclmodels


text:
      subroutine heapsort(face,n)
c
c  Ordina in senso crescente il vettore intero face
c
      implicit none
      integer  node,i,j,k,ik,jk
      integer  n
      integer  face(5,n)
      logical  conf_fac
c
      do node = 2,n
         i = node
         j = i/2
         do while((i.ne.1).and.(conf_fac(face(1,j),face(1,i))))
            call swap_t(face(1,j),face(1,i))
            i = j
            j = i/2
         end do
      end do
c
      do i = n,2,-1
         call swap_t(face(1,i),face(1,1))
         k = i-1
         ik = 1
         jk = 2
         if ((k.ge.3).and.(conf_fac(face(1,2),face(1,3))))  jk = 3
         do while ((jk.le.k).and.(conf_fac(face(1,ik),face(1,jk))))
            call swap_t(face(1,jk),face(1,ik))
            ik = jk
            jk = ik*2
            if ((jk+1.le.k).and.(conf_fac(face(1,jk),face(1,jk+1))))
     1         jk = jk+1
         end do
      end do
c
      return
c
      end


text:
      subroutine mapc2m_torus(blockno,xc,yc,xp,yp,zp,alpha)
      implicit none

      integer blockno
      double precision xc,yc,xp,yp,zp
      double precision alpha, r

      double precision pi

      common /compi/ pi

      r = 1 + alpha*cos(2*pi*yc)

      xp = r*cos(2*pi*xc)
      yp = r*sin(2*pi*xc)
      zp = alpha*sin(2*pi*yc)

      end

max_stars_repo_name: MelodyShih/forestclaw

max_issues_repo_name: scottaiton/forestclaw

max_forks_repo_name: scottaiton/forestclaw


text:
Schorre is an artist and jewelry designer in Davis, California.  Her pieces are created with natural and semiprecious stones, genuine Swarovski crystals, and sterling silver, gold filled and copper findings.  All pieces are oneofakind.  Prices range from $20 to $200.

Purchase  necklaces, bracelets & earrings  http://www.etsy.com/shop/oldcheviproductions Old Chevi Productions shopping cart.

Connect with http://www.facebook.com/pages/DavisCA/OldCheviProductions/151485727685 Old Chevi Productions on Facebook. http://twitter.com/OldChevi OldChevi is also using Twitter. http://instagram.com/schorreoldham OldChevi is also using Instagram.


text:
*DECK ALNREL
      FUNCTION ALNREL (X)
C***BEGIN PROLOGUE  ALNREL
C***PURPOSE  Evaluate ln(1+X) accurate in the sense of relative error.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C4B
C***TYPE      SINGLE PRECISION (ALNREL-S, DLNREL-D, CLNREL-C)
C***KEYWORDS  ELEMENTARY FUNCTIONS, FNLIB, LOGARITHM
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C ALNREL(X) evaluates ln(1+X) accurately in the sense of relative
C error when X is very small.  This routine must be used to
C maintain relative error accuracy whenever X is small and
C accurately known.
C
C Series for ALNR       on the interval -3.75000D-01 to  3.75000D-01
C                                        with weighted error   1.93E-17
C                                         log weighted error  16.72
C                               significant figures required  16.44
C                                    decimal places required  17.40
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CSEVL, INITS, R1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770401  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C***END PROLOGUE  ALNREL
      DIMENSION ALNRCS(23)
      LOGICAL FIRST
      SAVE ALNRCS, NLNREL, XMIN, FIRST
      DATA ALNRCS( 1) /   1.0378693562 743770E0 /
      DATA ALNRCS( 2) /   -.1336430150 4908918E0 /
      DATA ALNRCS( 3) /    .0194082491 35520563E0 /
      DATA ALNRCS( 4) /   -.0030107551 12753577E0 /
      DATA ALNRCS( 5) /    .0004869461 47971548E0 /
      DATA ALNRCS( 6) /   -.0000810548 81893175E0 /
      DATA ALNRCS( 7) /    .0000137788 47799559E0 /
      DATA ALNRCS( 8) /   -.0000023802 21089435E0 /
      DATA ALNRCS( 9) /    .0000004164 04162138E0 /
      DATA ALNRCS(10) /   -.0000000735 95828378E0 /
      DATA ALNRCS(11) /    .0000000131 17611876E0 /
      DATA ALNRCS(12) /   -.0000000023 54670931E0 /
      DATA ALNRCS(13) /    .0000000004 25227732E0 /
      DATA ALNRCS(14) /   -.0000000000 77190894E0 /
      DATA ALNRCS(15) /    .0000000000 14075746E0 /
      DATA ALNRCS(16) /   -.0000000000 02576907E0 /
      DATA ALNRCS(17) /    .0000000000 00473424E0 /
      DATA ALNRCS(18) /   -.0000000000 00087249E0 /
      DATA ALNRCS(19) /    .0000000000 00016124E0 /
      DATA ALNRCS(20) /   -.0000000000 00002987E0 /
      DATA ALNRCS(21) /    .0000000000 00000554E0 /
      DATA ALNRCS(22) /   -.0000000000 00000103E0 /
      DATA ALNRCS(23) /    .0000000000 00000019E0 /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  ALNREL
      IF (FIRST) THEN
         NLNREL = INITS (ALNRCS, 23, 0.1*R1MACH(3))
         XMIN = -1.0 + SQRT(R1MACH(4))
      ENDIF
      FIRST = .FALSE.
C
      IF (X .LE. (-1.0)) CALL XERMSG ('SLATEC', 'ALNREL', 'X IS LE -1',
     +   2, 2)
      IF (X .LT. XMIN) CALL XERMSG ('SLATEC', 'ALNREL',
     +   'ANSWER LT HALF PRECISION BECAUSE X TOO NEAR -1', 1, 1)
C
      IF (ABS(X).LE.0.375) ALNREL = X*(1. -
     1  X*CSEVL (X/.375, ALNRCS, NLNREL))
      IF (ABS(X).GT.0.375) ALNREL = LOG (1.0+X)
C
      RETURN
      END


! { dg-do run }
! PR 47065 - make sure that impure functions are not evaluated twice when
! replacing calls to trim with expression(1:len_trim)
module foo
  implicit none
contains
  function f()
    integer :: f
    integer :: s=0
    s = s + 1
    f = s
  end function f
end module foo

program main
  use foo
  implicit none
  character(len=10) :: line
  character(len=4) :: b(2)
  b(1) = 'a'
  b(2) = 'bc'
  write(unit=line,fmt='(A,A)') trim(b(f())), "X"
  if (line /= "aX          ") call abort
  if (f() .ne. 2) call abort
end program main


MODULE mod_mdct
IMPLICIT NONE
PRIVATE
PUBLIC :: mdct_initialize, sub_mdct
REAL (KIND = 8):: pi
INTEGER :: indx9(9)
COMPLEX (KIND = 8), SAVE:: omega9(0:8), omega9s(9), sqrt3_2, omega3(0:2), omega3s(3)
REAL (KIND = 8), PARAMETER :: c(0:7) = (/ -0.6000d0, -0.5350d0, -0.3300d0, -0.1850d0, -0.0950d0, -0.0410d0, -0.0142d0, -0.0037d0/)
REAL (KIND = 8), SAVE :: ca(0:7), cs(0:7), window_n(36), window_s(12), window_start(36), window_stop(36)
REAL (KIND = 8), SAVE, ALLOCATABLE :: subbuff(:, :, :)
CONTAINS
!---------------------------------------------------------------------------------------------
SUBROUTINE mdct_initialize()
IMPLICIT NONE
INTEGER :: i
pi = 4.0d0 * ATAN(1.0d0)
subbuff = 0.0d0
CALL fft9_initialize()
CALL fft3_initialize()
! ISO Table B.9 coefficients for aliasing reduction:
DO i = 0, 7
 cs(i) =  SQRT( 1.0d0    / ( 1.0d0 + c(i)**2 ) )
 ca(i) = -SQRT( c(i)**2  / ( 1.0d0 + c(i)**2 ) )
END DO
! ISO 2.4.3.4.10.3 Windowing, C.1.5.3.3
! normal window
DO i = 1, 36
 window_n(i) = SIN( REAL(2 * i - 1, KIND = 8) / 72.0d0 * pi)
END DO
! short window
DO i = 1, 12
 window_s(i) = SIN( REAL(2 * i - 1, KIND = 8) / 24.0d0 * pi)
END DO
! start / stop window
DO i = 1, 18
 window_start(i     ) = SIN( REAL(2 *  i       - 1, KIND = 8) / 72.0d0 * pi)
 window_stop (i + 18) = SIN( REAL(2 * (i + 18) - 1, KIND = 8) / 72.0d0 * pi)
END DO
DO i = 1, 6
 window_start(i + 24) = SIN( REAL(2 * (i + 6) - 1, KIND = 8) / 24.0d0 * pi)
 window_stop (i +  6) = SIN( REAL(2 *  i      - 1, KIND = 8) / 24.0d0 * pi)
END DO
window_start(19:24) = 1.0d0
window_start(31:36) = 0.0d0
window_stop ( 1: 6) = 0.0d0
window_stop (13:18) = 1.0d0
RETURN
END SUBROUTINE mdct_initialize
!-------------------------------------------------------------------------------------------------
SUBROUTINE cp_subband(subband, subbuff)
IMPLICIT NONE
REAL (KIND = 8), INTENT(IN ) :: subband(:, :, :) ! 32 * 36 * nchannel
REAL (KIND = 8), INTENT(OUT) :: subbuff(:, :, :) ! 32 * 54 * nchannel
subbuff = EOSHIFT(subbuff, 36, 0.0d0, 2)
subbuff(:     , 19:54  , :) =  subband(:, :, :)
subbuff(2:32:2, 20:54:2, :) = -subbuff(2:32:2, 20:54:2, :) ! ISO Figure A.4 Layer III decoder diagram
RETURN                                                     !    ~~~~~~~~~~~~ not written in the text but only in Figure
END SUBROUTINE cp_subband
!--------------------------------------------------------------------------------------------------
SUBROUTINE sub_mdct(subband, r_mdct, mblock_type, q_alias)
IMPLICIT NONE
REAL (KIND = 8), INTENT(IN ) :: subband(:, :, :)
REAL (KIND = 8), INTENT(OUT) :: r_mdct (:, :, :, :)
INTEGER        , INTENT(IN ) :: mblock_type(:)
LOGICAL        , INTENT(IN ) :: q_alias
INTEGER :: igranule, ichannel
LOGICAL, SAVE :: qfirst = .TRUE.
IF (qfirst) THEN
 qfirst = .FALSE.
 ALLOCATE( subbuff(32, 54, SIZE(subband, 3)) )
 subbuff = 0.0d0
 CALL mdct_initialize()
END IF
CALL cp_subband(subband, subbuff)
DO ichannel = 1, SIZE(subband, 3)
 DO igranule = 1, 2
  SELECT CASE (mblock_type(igranule))
   CASE ( 0) ! long block
    CALL sub_mdct_normal(igranule, 1, 32, subbuff(:, :, ichannel), r_mdct(:, :, igranule, ichannel))
    CALL anti_alias(r_mdct(1:32, :, igranule, ichannel))
   CASE (10) ! start block : long -> short
    CALL sub_mdct_start (igranule, 1, 32, subbuff(:, :, ichannel), r_mdct(:, :, igranule, ichannel))
    CALL anti_alias(r_mdct(1:32, :, igranule, ichannel))
   CASE (11) ! start block : long -> mixed
    CALL sub_mdct_normal(igranule, 1,  2, subbuff(:, :, ichannel), r_mdct(:, :, igranule, ichannel))
    CALL sub_mdct_start (igranule, 3, 32, subbuff(:, :, ichannel), r_mdct(:, :, igranule, ichannel))
    CALL anti_alias(r_mdct(1:32, :, igranule, ichannel))
   CASE (30) ! stop block : short -> long
    CALL sub_mdct_stop  (igranule, 1, 32, subbuff(:, :, ichannel), r_mdct(:, :, igranule, ichannel))
    CALL anti_alias(r_mdct(1:32, :, igranule, ichannel))
   CASE (31) ! stop block : mixed -> long
    CALL sub_mdct_normal(igranule, 1,  2, subbuff(:, :, ichannel), r_mdct(:, :, igranule, ichannel))
    CALL sub_mdct_stop  (igranule, 3, 32, subbuff(:, :, ichannel), r_mdct(:, :, igranule, ichannel))
    CALL anti_alias(r_mdct(1:32, :, igranule, ichannel))
   CASE (20) ! short block : short mode
    CALL sub_mdct_short(igranule, 1, 32, subbuff(:, :, ichannel), r_mdct(:, :, igranule, ichannel))
   CASE (21) ! short block : mixed mode
    CALL sub_mdct_normal(igranule, 1,  2, subbuff(:, :, ichannel), r_mdct(:, :, igranule, ichannel))
    CALL sub_mdct_short (igranule, 3, 32, subbuff(:, :, ichannel), r_mdct(:, :, igranule, ichannel))
    IF (q_alias) CALL anti_alias(r_mdct(1:2, :, igranule, ichannel))
   CASE DEFAULT
    STOP 'error : SUBROUTINE sub_mdct : block_type error'
  END SELECT
 END DO
END DO
RETURN
END SUBROUTINE sub_mdct
!--------------------------------------------------------------------------------------------------
SUBROUTINE sub_mdct_normal(igranule, n0, n1, subbuff, r_mdct)
IMPLICIT NONE
INTEGER        , INTENT(IN ) :: igranule, n0, n1
REAL (KIND = 8), INTENT(IN ) :: subbuff(:, :)
REAL (KIND = 8), INTENT(OUT) :: r_mdct (:, :)
REAL (KIND = 8)              :: wk(36)
INTEGER :: iband, k
k = 18 * (igranule - 1) + 1
DO iband = n0, n1
 wk = window_n * subbuff(iband, k:k + 35)
 CALL mdct36(wk, r_mdct(iband, :))
END DO
RETURN
END SUBROUTINE sub_mdct_normal
!--------------------------------------------------------------------------------------------------
SUBROUTINE sub_mdct_start(igranule, n0, n1, subbuff, r_mdct)
IMPLICIT NONE
INTEGER        , INTENT(IN ) :: igranule, n0, n1
REAL (KIND = 8), INTENT(IN ) :: subbuff(:, :)
REAL (KIND = 8), INTENT(OUT) :: r_mdct (:, :)
REAL (KIND = 8)              :: wk(36)
INTEGER :: iband, k
k = 18 * (igranule - 1) + 1
DO iband = n0, n1
 wk = window_start * subbuff(iband, k:k + 35)
 CALL mdct36(wk, r_mdct(iband, :))
END DO
RETURN
END SUBROUTINE sub_mdct_start
!--------------------------------------------------------------------------------------------------
SUBROUTINE sub_mdct_stop(igranule, n0, n1, subbuff, r_mdct)
IMPLICIT NONE
INTEGER        , INTENT(IN ) :: igranule, n0, n1
REAL (KIND = 8), INTENT(IN ) :: subbuff(:, :)
REAL (KIND = 8), INTENT(OUT) :: r_mdct (:, :)
REAL (KIND = 8)              :: wk(36)
INTEGER :: iband, k
k = 18 * (igranule - 1) + 1
DO iband = n0, n1
 wk = window_stop * subbuff(iband, k:k + 35)
 CALL mdct36(wk, r_mdct(iband, :))
END DO
RETURN
END SUBROUTINE sub_mdct_stop
!--------------------------------------------------------------------------------------------------
SUBROUTINE sub_mdct_short(igranule, n0, n1, subbuff, r_mdct)
IMPLICIT NONE
INTEGER        , INTENT(IN ) :: igranule, n0, n1
REAL (KIND = 8), INTENT(IN ) :: subbuff(:, :)
REAL (KIND = 8), INTENT(OUT) :: r_mdct (:, :)
REAL (KIND = 8)              :: wk(12)
INTEGER :: iband, k, m, m0, n
m0 = 18 * (igranule - 1) + 1 + 6
DO k = 1, 3
 m = m0 + 6 * (k - 1)
 n =  1 + 6 * (k - 1)
 DO iband = n0, n1
  wk = window_s * subbuff(iband, m:m + 11)
  CALL mdct12(wk, r_mdct(iband, n:n + 5))
 END DO
END DO
RETURN
END SUBROUTINE sub_mdct_short
!--------------------------------------------------------------------------------------------------
SUBROUTINE anti_alias(x) ! C.1.5.3.3 Aliasing-Butterfly, 2.4.3.4.10.1 Alias Reduction, Table B.9
IMPLICIT NONE
REAL (KIND = 8), INTENT(IN OUT) :: x(:, :) ! x(32, 18)
REAL (KIND = 8) :: tmp1, tmp2
INTEGER :: iband, k
DO iband = 1, SIZE(x, 1) - 1
 DO k = 0, 7
  tmp1 = x(iband    , 18 - k)
  tmp2 = x(iband + 1,  1 + k)
  x(iband    , 18 - k) =  tmp1 * cs(k) + tmp2 * ca(k)
  x(iband + 1,  1 + k) =  tmp2 * cs(k) - tmp1 * ca(k)
 END DO
END DO
RETURN
END SUBROUTINE anti_alias
!--------------------------------------------------------------------------------------------------
SUBROUTINE mdct36(x_in, x_out) ! ISO C.1.5.3.3 MDCT:   ! MDCT(4n) <=> FFT(n) proof: http://members.at.infoseek.co.jp/kitaurawa/mdct.pdf   (JP)
IMPLICIT NONE                                          !                            http://members.at.infoseek.co.jp/kitaurawa/mdct_e.pdf (E)
REAL (KIND = 8), INTENT(IN ) :: x_in (:) ! 36
REAL (KIND = 8), INTENT(OUT) :: x_out(:) ! 18
REAL (KIND = 8):: x_shift(36), x_re, x_im
COMPLEX (KIND = 8):: fft(9)
INTEGER:: i
x_shift( 1: 9) = -x_in(28:36)
x_shift(10:36) =  x_in( 1:27)
DO i = 1, 9
 x_re = x_shift(2 * i -  1) - x_shift(38 - 2 * i)
 x_im = x_shift(2 * i + 17) - x_shift(20 - 2 * i)
 fft(i) = CMPLX(x_re, x_im, KIND = 8) * omega9s(i)
END DO
CALL fft9(fft)
fft = fft * omega9s
DO i = 1, 9
 x_out(2 * i - 1) =  REAL(      fft(i)      , KIND = 8)
 x_out(2 * i    ) =  REAL(AIMAG(fft(10 - i)), KIND = 8)
END DO
RETURN
END SUBROUTINE mdct36
!---------------------------------------------------------------------------------------------------
SUBROUTINE fft9(fft) ! FFT of 3^n (n=2)
IMPLICIT NONE
COMPLEX (KIND = 8), INTENT(IN OUT):: fft(:)
COMPLEX (KIND = 8) :: c1, c2, c3, c4, tmp1, tmp2, tmp3
INTEGER, PARAMETER :: nn = 9
INTEGER :: i, j, iphase1, iphase2, m1, m2, m3, k3, kn3
fft = fft(indx9) / 9.0d0 ! reorder and normalize
DO k3 = 1, 2 ! 3^n (n=2)
 kn3 = 3**(k3 - 1)
 DO i = 1, nn, 3 * kn3
  DO j = 1, kn3
   iphase1 = 3**(2 - k3) * (j - 1)
   iphase2 = 2 * iphase1
   c1 = omega9( MOD(iphase1, nn) )
   c2 = omega9( MOD(iphase2, nn) )
   m1 = i + j - 1
   m2 = m1 + kn3
   m3 = m2 + kn3
   tmp1 =      fft(m1)
   tmp2 = c1 * fft(m2)
   tmp3 = c2 * fft(m3)
   fft(m1) = tmp1 + tmp2 + tmp3
   c3 = tmp1 - 0.5d0 * ( tmp2 + tmp3 )
   c4 =     sqrt3_2  * ( tmp2 - tmp3 ) ! sqrt3_2 = i sqrt(3) / 2
   fft(m2) = c3 + c4
   fft(m3) = c3 - c4
  END DO
 END DO
END DO
RETURN
END SUBROUTINE fft9
!-----------------------------------------------------------------------
SUBROUTINE fft9_initialize()
IMPLICIT NONE
INTEGER :: i, j, k, n
pi = 4.0d0 * ATAN(1.0d0)
sqrt3_2 = CMPLX(0.0d0, SQRT(0.75d0), KIND = 8) ! iSQRT(3) / 2
DO i = 1, 9
 n = 0
 k = i - 1
 DO j = 2 - 1, 0, -1
  n = n + MOD(k, 3) * 3**j
  k = k / 3
 END DO
 indx9(i) = n + 1
 !omega9(i - 1) = EXP( CMPLX(0.0d0, 2.0d0 * pi /  9.0d0 * REAL(i - 1, KIND = 8), KIND = 8) )
 !omega9s(i)    = EXP( CMPLX(0.0d0, 2.0d0 * pi / 36.0d0 * ( 1.0d0 / 8.0d0 + REAL(i - 1, KIND = 8) ), KIND = 8) )
 omega9(i - 1) = CMPLX( COS( REAL( 2 * i - 2, KIND = 8) /   9.0d0 * pi ),   &
                        SIN( REAL( 2 * i - 2, KIND = 8) /   9.0d0 * pi ), KIND = 8 )
 omega9s(i)    = CMPLX( COS( REAL( 8 * i - 7, KIND = 8) / 144.0d0 * pi ), &
                        SIN( REAL( 8 * i - 7, KIND = 8) / 144.0d0 * pi ), KIND = 8 )
END DO
RETURN
END SUBROUTINE fft9_initialize
!-----------------------------------------------------------------------
SUBROUTINE mdct12(x_in, x_out)
IMPLICIT NONE
REAL (KIND = 8), INTENT(IN ) :: x_in (:) ! 12
REAL (KIND = 8), INTENT(OUT) :: x_out(:) !  6
REAL (KIND = 8):: x_shift(12), x_re, x_im
COMPLEX (KIND = 8):: fft(3)
INTEGER:: i
x_shift( 1: 3) = -x_in(10:12)
x_shift( 4:12) =  x_in( 1: 9)
DO i = 1, 3
 x_re = x_shift(2 * i - 1) - x_shift(14 - 2 * i)
 x_im = x_shift(2 * i + 5) - x_shift( 8 - 2 * i)
 fft(i) = CMPLX(x_re, x_im, KIND = 8) * omega3s(i)
END DO
CALL fft3(fft)
fft = fft * omega3s
DO i = 1, 3
 x_out(2 * i - 1) =  REAL(      fft(i)     , KIND = 8)
 x_out(2 * i    ) =  REAL(AIMAG(fft(4 - i)), KIND = 8)
END DO
RETURN
END SUBROUTINE mdct12
!-----------------------------------------------------------------------
SUBROUTINE fft3(fft)
IMPLICIT NONE
COMPLEX (KIND = 8), INTENT(IN OUT):: fft(:)
COMPLEX (KIND = 8) :: c2, c3, c4
c2 = fft(2) + fft(3)
c3 = fft(1) - 0.5d0 * c2
c4 = sqrt3_2 * ( fft(2) - fft(3) ) ! sqrt3_2 = i sqrt(3) / 2
fft(1) = fft(1) + c2
fft(2) = c3 + c4
fft(3) = c3 - c4
fft = fft / 3.0d0 ! normalization
RETURN
END SUBROUTINE fft3
!--------------------------------------------------------------------------------------------------
SUBROUTINE fft3_initialize()
IMPLICIT NONE
INTEGER :: i
pi = 4.0d0 * ATAN(1.0d0)
sqrt3_2 = CMPLX(0.0d0, SQRT(0.75d0), KIND = 8) ! iSQRT(3) / 2
DO i = 1, 3
 !omega3(i - 1) = EXP( CMPLX(0.0d0, 2.0d0 * pi /  3.0d0 * REAL(i - 1, KIND = 8), KIND = 8) )
 !omega3s(i)    = EXP( CMPLX(0.0d0, 2.0d0 * pi / 12.0d0 * ( 1.0d0 / 8.0d0 + REAL(i - 1, KIND = 8) ), KIND = 8) )
 omega3(i - 1) = CMPLX( COS( REAL( 2 * i - 2, KIND = 8) /  3.0d0 * pi ),   &
                        SIN( REAL( 2 * i - 2, KIND = 8) /  3.0d0 * pi ), KIND = 8 )
 omega3s(i)    = CMPLX( COS( REAL( 8 * i - 7, KIND = 8) / 48.0d0 * pi ),   &
                        SIN( REAL( 8 * i - 7, KIND = 8) / 48.0d0 * pi ), KIND = 8 )
END DO
RETURN
END SUBROUTINE fft3_initialize
!-----------------------------------------------------------------------
END MODULE mod_mdct


text:
c
c   Dint – version 2.0  is licensed under the Apache License, Version 2.0 (the "License");
c   you may not use Dint – version 2.0 except in compliance with the License.
c   You may obtain a copy of the License at
c       http://www.apache.org/licenses/LICENSE-2.0
c   The license is also given in the LICENSE file.
c   Unless required by applicable law or agreed to in writing, software
c   distributed under the License is distributed on an "AS IS" BASIS,
c   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
c   See the License for the specific language governing permissions and limitations under the License.
c
c -------------------------------------------------------------------------------------------
c  Dint : Direct Nonadiabatic Trajectories A code for non-Born–Oppenheimer molecular dynamics
c
c  version 2.0
c
c  A. W. Jasper
c  Argonne National Laboratory
c
c  Rui Ming Zhang
c  Tsinghua University
c
c  and
c
c  D. G. Truhlar
c  University of Minnesota
c
c  copyright  2020
c  Donald G. Truhlar and Regents of the University of Minnesota
c----------------------------------------------------------------------------------------------

      subroutine header

      write(6,*)" **************************************************"
      write(6,*)" **                     Dint                     **"
      write(6,*)" **       Direct Nonadiabatic Trajectories       **"
      write(6,*)" **                                              **"
      write(6,*)" **                 Version 2.0                  **"
      write(6,*)" **              December 16, 2020               **"
      write(6,*)" **                                              **"
      write(6,*)" **                A. W. Jasper                  **"
      write(6,*)" **         Sandia National Laboratories         **"
      write(6,*)" **               Rui Ming Zhang                 **"
      write(6,*)" **            Tsinghua University               **"
      write(6,*)" **                     and                      **"
      write(6,*)" **                D. G. Truhlar                 **"
      write(6,*)" **           University of Minnesota            **"
      write(6,*)" **************************************************"
      write(6,*)

      return
      end

max_stars_repo_name: kevinmooreiii/DiNT

max_issues_repo_name: kevinmooreiii/DiNT

max_forks_repo_name: kevinmooreiii/DiNT


text:
!********************************************************
!
!     wrapNNPDF.f:
!     Routine called by LHAPDF package for the evolved PDF
!     in a (x,Q) point as called by NNPDF.LHpdf file.
!
!     In 'wrapevolve.f' the package calls:
!     IF(NAME(NSET).EQ.'NNPDF') call NNPDFevolve(x,Q,f)
!     IF(NAME(NSET).EQ.'NNPDF') call NNPDFread(nset)
!     IF(NAME(NSET).EQ.'NNPDF') call NNPDFalfa(alfas,Q)
!     IF(NAME(NSET).EQ.'NNPDF') call NNPDFinit(nset,Eorder,Q2fit)
!
!
!********************************************************

      SUBROUTINE NNPDFevolve(x,Q,f)
      IMPLICIT none
!
      INCLUDE 'parmsetup.inc'
      CHARACTER*16 name(nmxset)
      INTEGER nmem(nmxset),ndef(nmxset),mmem
      COMMON/NAME/name,nmem,ndef,mmem
!
      INTEGER nset,pdfmem
      REAL*8 f(-6:6)
      REAL*8 x,Q
!
      INTEGER order
      REAL*8 alfas,alphaNNPDF,alfas0,q0
      REAL*8 Eorder,Q2fit,mass
      REAL*8 pdfout(13)
      REAL*8 check
      REAL*8 pi
      PARAMETER (pi = 3.1415926535897932385)
!
      INTEGER i,nff
      REAL*8 xpt
!
      INTEGER ipt,imodev,ivfn,itmc
      COMMON/nnpdf10EVFLAGS/ipt,imodev,ivfn,itmc
      INTEGER ieval,niter,nmax
      COMMON/nnpdf10ngrid/ ieval,niter,nmax
      REAL*8 xmin,xm1,xm2,xmax
      COMMON/nnpdf10GRID/xmin,xm1,xm2,xmax
      REAL*8 epsevol
      COMMON/nnpdf10evolacc/epsevol
!
      REAL*8 qq,qq2,qth(4:6)
      REAL*8 q20,q2
      COMMON/nnpdf10EVSCALE/q20,q2
      REAL*8 q2th(4:6),asref,q2ref
      COMMON/nnpdf10vfns/q2th,asref,q2ref
      REAL*8 as0,asc,asb,ast,asq
      COMMON/nnpdf10AS/ as0,asc,asb,ast,asq
!
      q2    = q**2.                ! for the commons
      qq2   = q**2.
      qq = q
      asq   = alphaNNPDF(qq)
      asq   = asq / (4.*pi)
!
      xpt = x
      CALL lh_evolfactx(0,xpt)
      CALL lh_pdfevolx(xpt,pdfout)
      CALL lh_pdfevln2lha(pdfout,f)
      DO i = -6,6,1
         f(i) = x * f(i)
      ENDDO
!
      RETURN

!********************************************************

      ENTRY NNPDFread(nset)
!
!     Read evol params from NNPDF08.LHpdf
      READ(1,*) ivfn,imodev,itmc,niter,xmin,xm1,xm2,xmax,epsevol
!
      RETURN

!********************************************************

      ENTRY NNPDFalfa(alfas,Q)
!
      QQ = Q
      alfas =  alphaNNPDF(QQ)
!
      RETURN

!********************************************************

      ENTRY NNPDFinit(nset,Eorder,Q2fit)
!
      CALL GetOrderPDFM(nset,order)
      IPT = order
!
      CALL GetQ2fitM(nset,Q2)
      Q2fit = Q2
      Q20   = Q2
!
      call GetQmassM(nset,4,mass)
      QTH(4) = mass
      call GetQmassM(nset,5,mass)
      QTH(5) = mass
      call GetQmassM(nset,6,mass)
      QTH(6) = mass
      DO i = 4,6
         q2th(i) = qth(i)**2d0
      ENDDO
!
      CALL GetAlfas(nset,alfas0,Q0)
      asref = 0.119
      q2ref = q0**2d0
!
      RETURN

!********************************************************

      ENTRY NNPDFpdf(nset)
!
      pdfmem = mmem
      IF (pdfmem.LT.0) THEN
         WRITE(*,*) 'NNPDF set:'
         WRITE(*,*) 'PDF member out of range:'
         WRITE(*,*) 'member = ',pdfmem
         STOP
      ENDIF
      RETURN
!
      END

max_stars_repo_name: AllaMaevskaya/AliRoot

max_issues_repo_name: AllaMaevskaya/AliRoot

max_forks_repo_name: AllaMaevskaya/AliRoot