!-----------------------------------------------------------------------
!
! RAdS Jan12.
! Copyright (C) 2011-2014 by Rodrigo Alvares de Souza.
! Mail: <rsouza01@gmail.com>. Web: "http://www.astro.iag.usp.br/~rsouza/".
! This program may be copied and/or distributed freely. See the
! _ terms and conditions in the files in the doc/ subdirectory.
!
!
! MODULE: tov_module
!
!> @author
!> Rodrigo Souza
!
! DESCRIPTION:
!> TOV related types and functions (EDO's and stuff)
!
!-----------------------------------------------------------------------
module tov_module

use types
use global_constants
use config_module
use output_module
use cubic_spline_module
use stat_module

    implicit none

    contains

    !> \brief Set initial conditions
    !!
    !! \param parameters
    !! \param Y
    !! \param V
    !!
    subroutine set_initial_conditions(parameters, Y, V)
    implicit none

        type(ConfigParameters) :: parameters
        double precision :: Y(N_DIFF_EQUATIONS)
        double precision :: V(N_VARIABLES)

        !Initial conditions

        !n = R_bar = 0
        V(IDX_RADIUS_BAR) = 0.

        !M_bar(0) = 0
        V(IDX_MASS_BAR) = 0.

        !Rho_bar(0) = RHO_0 / RHO_ADIM
        V(IDX_RHO_BAR) = parameters%RHO_0 / parameters%RHO_ADIM

        !P_bar(0) = P_0/(Rho_0 * c^2)
        V(IDX_PRESSURE_BAR) = parameters%P_0_bar

        V(IDX_ENERGY_DENSITY_BAR) = V(IDX_RHO_BAR)

        V(IDX_BARYON_NUMBER) = 0.

        V(IDX_INFOR_ENTROPY) = 0.
        V(IDX_DISEQUILIBRIUM) = 0.

        Y(IDX_EQN_DM_DN) = 0.
        Y(IDX_EQN_DP_DN) = V(IDX_PRESSURE_BAR)

        !write (*,*) &
        !    '(IDX_MASS_BAR, IDX_PRESSURE_BAR, IDX_ENERGY_DENSITY_BAR, IDX_RADIUS_BAR, IDX_RHO_BAR) = (', &
        !    V(IDX_MASS_BAR), ' ,', V(IDX_PRESSURE_BAR), ' ,', V(IDX_ENERGY_DENSITY_BAR), ' ,', &
        !    V(IDX_RADIUS_BAR), ' ,', V(IDX_RHO_BAR), ')'

    end subroutine set_initial_conditions

    !> \brief The calculations that should be performed after the runge-kutta.
    !! Conversions and other values that should be calculated go here
    !!
    !! \param parameters
    !! \param t
    !! \param Y
    !! \param V
    !!
    subroutine process_specific_calculations(parameters, t, Y, V)
        type(ConfigParameters) :: parameters
        double precision :: t
        double precision :: Y(N_DIFF_EQUATIONS)
        double precision :: V(N_VARIABLES)

        !OK
        !write (*,*) &
        !    '(IDX_MASS_BAR, IDX_PRESSURE_BAR, IDX_ENERGY_DENSITY_BAR, IDX_RADIUS_BAR, IDX_RHO_BAR) = (', &
        !    V(IDX_MASS_BAR), ' ,', V(IDX_PRESSURE_BAR), ' ,', V(IDX_ENERGY_DENSITY_BAR), ' ,', &
        !    V(IDX_RADIUS_BAR), ' ,', V(IDX_RHO_BAR), ')'

        V(IDX_RADIUS_BAR) = t;

        V(IDX_PRESSURE_BAR) = Y(IDX_EQN_DP_DN);

        V(IDX_ENERGY_DENSITY_BAR) = energy_density_from_eos(parameters, t, Y, V);

        V(IDX_RHO_BAR) = V(IDX_ENERGY_DENSITY_BAR);

        V(IDX_MASS_BAR) = Y(IDX_EQN_DM_DN);

        V(IDX_INFOR_ENTROPY) = V(IDX_INFOR_ENTROPY) + calc_infor_entropy(parameters, t, Y, V);
        V(IDX_DISEQUILIBRIUM) = V(IDX_DISEQUILIBRIUM) + calc_infor_disequilibrium(parameters, t, Y, V);

        V(IDX_SPEED_OF_SOUND) = SQRT((parameters%P_0 * V(IDX_PRESSURE_BAR))/ &
            (parameters%ENERGY_DENSITY_0 * V(IDX_RHO_BAR))) * &
            LIGHT_SPEED_SCALE(parameters%UNIT_SYSTEM)

        if (parameters%eos_file_provides_baryonic_density) then
            V(IDX_BARYON_NUMBER) = V(IDX_BARYON_NUMBER) + calc_baryon_number(parameters, t, Y, V);
        end if

        !write (*,*) V(IDX_ENERGY_DENSITY_BAR) * (parameters%RHO_ADIM * LIGHT_SPEED(parameters%UNIT_SYSTEM)**2.), &
        !', ', V(IDX_PRESSURE_BAR)* (parameters%RHO_ADIM * LIGHT_SPEED(parameters%UNIT_SYSTEM)**2.)

    end subroutine process_specific_calculations

    !---------------------------------------------------------------------------
    ! DESCRIPTION:
    !> Evaluates derivatives at a point
    !---------------------------------------------------------------------------
    double precision function derivative(parameters, t, Y, V, equationNumber)
    implicit none
        type(ConfigParameters) :: parameters
        double precision :: t, Y(N_DIFF_EQUATIONS), V(N_VARIABLES)
        integer :: equationNumber

        double precision :: returnValue

        !write (*,*) &
        !    '(IDX_MASS_BAR, IDX_PRESSURE_BAR, IDX_ENERGY_DENSITY_BAR, IDX_RADIUS_BAR, IDX_RHO_BAR) = (', &
        !    V(IDX_MASS_BAR), ' ,', V(IDX_PRESSURE_BAR), ' ,', V(IDX_ENERGY_DENSITY_BAR), ' ,', &
        !    V(IDX_RADIUS_BAR), ' ,', V(IDX_RHO_BAR), ')'

        select case (equationNumber)
            case (IDX_EQN_DM_DN)
                  returnValue = dM_dn(parameters, t, Y, V);

            case (IDX_EQN_DP_DN)

                    if (parameters%RHO_0 <= parameters%cutoff_RHO_0) then
                        returnValue = dP_dN_Newtonian(parameters, t, Y, V);
                    else
                        returnValue = dP_dn_Relativistic(parameters, t, Y, V);
                    end if

            case default
                returnValue = 0.;

        end select

        derivative = returnValue

    end function derivative

    !---------------------------------------------------------------------------
    ! DESCRIPTION:
    !> Evaluates the pressure P at a point
    !
    !> Compute \f$ P \f$
    !---------------------------------------------------------------------------
    double precision function dM_dn(parameters, t, Y, V)
    implicit none
        type(ConfigParameters) :: parameters
        double precision :: t
        double precision :: Y(N_DIFF_EQUATIONS)
        double precision :: V(N_VARIABLES)

        double precision :: energy_density_bar

        !write (*,*) '(t, energy density) = (', t, ', ', V(IDX_ENERGY_DENSITY_BAR), ')'

        !V(IDX_ENERGY_DENSITY_BAR) = energy_density_from_eos(parameters, t, Y, V)
        energy_density_bar = energy_density_from_eos(parameters, t, Y, V);

        dM_dn = energy_density_bar * t**2.

    end function dM_dn

    !---------------------------------------------------------------------------
    ! DESCRIPTION:
    !> Evaluates the pressure P at a point (TOV)
    !
    !> Compute \f$ P \f$
    !---------------------------------------------------------------------------
    double precision function dP_dn_Relativistic(parameters, t, Y, V)
    implicit none
        type(ConfigParameters) :: parameters
        double precision :: t
        double precision :: Y(N_DIFF_EQUATIONS)
        double precision :: V(N_VARIABLES)

        double precision :: numerator, denominator

        double precision :: energy_density_bar
        double precision :: pressure_bar

        energy_density_bar = energy_density_from_eos(parameters, t, Y, V);
        pressure_bar = Y(IDX_PRESSURE_BAR);

        !write (*,*) &
        !    '(t, MASS, PRESSURE, ENERGY, RADIUS, RHO) = (', t, ', ', &
        !    V(IDX_MASS_BAR), ' ,', V(IDX_PRESSURE_BAR), ' ,', V(IDX_ENERGY_DENSITY_BAR), ' ,', &
        !    V(IDX_RADIUS_BAR), ' ,', V(IDX_RHO_BAR), ')'

        numerator = - (energy_density_bar + pressure_bar)*(pressure_bar*(t**3.) + V(IDX_MASS_BAR));

        denominator = (t**2.) * (1. - 2. * V(IDX_MASS_BAR) / t);

        !write (*,*) '(t, numerator, denominator)', t, ', ', numerator, ', ', denominator, ')'

        dP_dn_Relativistic = numerator / denominator;
    end function dP_dn_Relativistic

    !---------------------------------------------------------------------------
    ! DESCRIPTION:
    !> Evaluates the pressure P at a point (Newtonian)
    !
    !> Compute \f$ P \f$
    !---------------------------------------------------------------------------
    double precision function dP_dn_Newtonian(parameters, t, Y, V)
    implicit none
        type(ConfigParameters) :: parameters
        double precision :: t, Y(N_DIFF_EQUATIONS), V(N_VARIABLES)

        dP_dn_Newtonian = - (V(IDX_MASS_BAR) * V(IDX_ENERGY_DENSITY_BAR) / (t **2.))
    end function dP_dn_Newtonian

    !> \brief Function that decides whether the program must stop.
    !!
    !! \param parameters
    !! \param t
    !! \param Y
    !! \param V
    !! \return
    !!
    logical function can_stop(parameters, t, Y, V)
    implicit none
        type(ConfigParameters) :: parameters
        double precision :: t
        double precision :: Y(N_DIFF_EQUATIONS)
        double precision :: V(N_VARIABLES)

        !can_stop = (V(IDX_PRESSURE_BAR) <= 0.0 .OR. V(IDX_ENERGY_DENSITY_BAR) <= 0.0);
        can_stop = (V(IDX_RHO_BAR) <= parameters%cutoff_density_bar .OR. &
            V(IDX_PRESSURE_BAR) <= 0.0);

        !DEBUG
        !write (*,*) 'cutoff_density_bar, V(IDX_RHO_BAR) => ', parameters%cutoff_density_bar, &
        !    ',', V(IDX_RHO_BAR)

    end function can_stop

end module tov_module


text:
C
      SUBROUTINE FELDG(GLAT,GLON,ALT,BNORTH,BEAST,BDOWN,BABS)
c  ------------------------------------------------------------------------
c  Part of "poa_calfos" (STECF/IPMG)
c  Part of "NEWBILCAL" program to produce IGRF geomagnetic model (GSFC)
c
c  History:
c  --------
c  Version   Date        Author          Description
c      1.0       Jul 00  M. Rosa         POA calibration code
c  ------------------------------------------------------------------------
C-------------------------------------------------------------------
C CALCULATES EARTH MAGNETIC FIELD FROM SPHERICAL HARMONICS MODEL
C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTRE, INTERNAL NOTE 61,
C      1970.
C--------------------------------------------------------------------
C CHANGES (D. BILITZA, NOV 87):
C   - FIELD COEFFICIENTS IN BINARY DATA FILES INSTEAD OF BLOCK DATA
C   - CALCULATES DIPOL MOMENT
C--------------------------------------------------------------------
C  INPUT:  ENTRY POINT FELDG
C               GLAT  GEODETIC LATITUDE IN DEGREES (NORTH)
C               GLON  GEODETIC LONGITUDE IN DEGREES (EAST)
C               ALT   ALTITUDE IN KM ABOVE SEA LEVEL
C
C          ENTRY POINT FELDC
C              V(3)  CARTESIAN COORDINATES IN EARTH RADII (6371.2 KM)
C                     X-AXIS POINTING TO EQUATOR AT 0 LONGITUDE
C                     Y-AXIS POINTING TO EQUATOR AT 90 LONG.
C                     Z-AXIS POINTING TO NORTH POLE
C
C          COMMON BLANK AND ENTRY POINT FELDI ARE NEEDED WHEN USED
C            IN CONNECTION WITH L-CALCULATION PROGRAM SHELLG.
C
C          COMMON /MODEL/ AND /GENER/
C              UMR     = ATAN(1.0)*4./180.   <DEGREE>*UMR=<RADIANT>
C              ERA       EARTH RADIUS FOR NORMALIZATION OF CARTESIAN
C                     COORDINATES (6371.2 KM)
C              AQUAD, BQUAD   SQUARE OF MAJOR AND MINOR HALF AXIS FOR
C                     EARTH ELLIPSOID AS RECOMMENDED BY INTERNATIONAL
C                     ASTRONOMICAL UNION (6378.160, 6356.775 KM).
c comment (MRR): This is IAU 1964 (r=6378.160; 1/f=298.25)
C              NMAX    MAXIMUM ORDER OF SPHERICAL HARMONICS
C              TIME       YEAR (DECIMAL: 1973.5) FOR WHICH MAGNETIC
C                     FIELD IS TO BE CALCULATED
C              G(M)       NORMALIZED FIELD COEFFICIENTS (SEE FELDCOF)
C                     M=NMAX*(NMAX+2)
C------------------------------------------------------------------------
C  OUTPUT: BABS   MAGNETIC FIELD STRENGTH IN GAUSS
C          BNORTH, BEAST, BDOWN   COMPONENTS OF THE FIELD WITH RESPECT
C                TO THE LOCAL GEODETIC COORDINATE SYSTEM, WITH AXIS
C                POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
C                AND DOWNWARD.
C-----------------------------------------------------------------------

      INTEGER  NMAX,I,K,IH,IL,IS,IHMAX,LAST,M,IMAX
      REAL*4   V(3),B(3)
      REAL*4   TIME
      REAL*4   G(144)
      REAL*4   f,x,y,z,s,t,bnorth,beast,bdown,babs
      REAL*4   glat,glon,alt,rlat,ct,st,d,rlon,cp,sp,rq
      real*4   bxxx,byyy,bzzz,brho,xxx,yyy,zzz,rho
      REAL*4   UMR,ERAD,AQUAD,BQUAD
      REAL*4   XI(3),H(144)
      CHARACTER*33        NAME

      COMMON /SHCOM/  XI,H
c ** MR ** changed Common to align on boundaries
c ** org **     COMMON/MODEL/       NAME,NMAX,TIME,G
      COMMON/MODEL/       TIME,G,NMAX,NAME
      COMMON/GENER/       UMR,ERAD,AQUAD,BQUAD
C
C-- IS RECORDS ENTRY POINT
C
C*****ENTRY POINT  FELDG  TO BE USED WITH GEODETIC CO-ORDINATES
      IS=1
      RLAT=GLAT*UMR
      CT=SIN(RLAT)
      ST=COS(RLAT)
      D=SQRT(AQUAD-(AQUAD-BQUAD)*CT*CT)
      RLON=GLON*UMR
      CP=COS(RLON)
      SP=SIN(RLON)
      ZZZ=(ALT+BQUAD/D)*CT/ERAD
      RHO=(ALT+AQUAD/D)*ST/ERAD
      XXX=RHO*CP
      YYY=RHO*SP
      GOTO10
      ENTRY FELDC(V,B)
C*****ENTRY POINT  FELDC  TO BE USED WITH CARTESIAN CO-ORDINATES
      IS=2
      XXX=V(1)
      YYY=V(2)
      ZZZ=V(3)
10    RQ=1./(XXX*XXX+YYY*YYY+ZZZ*ZZZ)
      XI(1)=XXX*RQ
      XI(2)=YYY*RQ
      XI(3)=ZZZ*RQ
      GOTO20
      ENTRY FELDI
C*****ENTRY POINT  FELDI  USED FOR L COMPUTATION
      IS=3
20    IHMAX=NMAX*NMAX+1
      LAST=IHMAX+NMAX+NMAX
      IMAX=NMAX+NMAX-1
      DO 8 I=IHMAX,LAST
8     H(I)=G(I)
      DO 6 K=1,3,2
      I=IMAX
      IH=IHMAX
1     IL=IH-I
      F=2./FLOAT(I-K+2)
      X=XI(1)*F
      Y=XI(2)*F
      Z=XI(3)*(F+F)
      I=I-2
      IF(I-1)5,4,2
2     DO 3 M=3,I,2
      H(IL+M+1)=G(IL+M+1)+Z*H(IH+M+1)+X*(H(IH+M+3)-H(IH+M-1))
     &                               -Y*(H(IH+M+2)+H(IH+M-2))
3     H(IL+M)=G(IL+M)+Z*H(IH+M)+X*(H(IH+M+2)-H(IH+M-2))
     &                         +Y*(H(IH+M+3)+H(IH+M-1))
4     H(IL+2)=G(IL+2)+Z*H(IH+2)+X*H(IH+4)-Y*(H(IH+3)+H(IH))
      H(IL+1)=G(IL+1)+Z*H(IH+1)+Y*H(IH+4)+X*(H(IH+3)-H(IH))
5     H(IL)=G(IL)+Z*H(IH)+2.*(X*H(IH+1)+Y*H(IH+2))
      IH=IL
      IF(I.GE.K)GOTO1
6     CONTINUE
      IF(IS.EQ.3)RETURN
      S=.5*H(1)+2.*(H(2)*XI(3)+H(3)*XI(1)+H(4)*XI(2))
      T=(RQ+RQ)*SQRT(RQ)
      BXXX=T*(H(3)-S*XXX)
      BYYY=T*(H(4)-S*YYY)
      BZZZ=T*(H(2)-S*ZZZ)
      IF(IS.EQ.2)GOTO7
      BABS=SQRT(BXXX*BXXX+BYYY*BYYY+BZZZ*BZZZ)
      BEAST=BYYY*CP-BXXX*SP
      BRHO=BYYY*SP+BXXX*CP
      BNORTH=BZZZ*ST-BRHO*CT
      BDOWN=-BZZZ*CT-BRHO*ST
      RETURN
7     B(1)=BXXX
      B(2)=BYYY
      B(3)=BZZZ
      RETURN
      END
C

ext: f

max_stars_repo_name: iraf-community/stsdas


text:
! Copyright (C) 2011-2013 German Research School for Simulation Sciences GmbH,
!              Aachen and others.
! Please see the ../COPYRIGHT file one directory above for details.

!> A module providing quadruple number input to Lua functions
!!
!! Note that Lua actually only handles double precision, and the numbers are
!! converted accordingly. Thus this is merely a convenience interface, to allow
!! the usage of the functions from this module with quadruple precision numbers.
module aot_quadruple_fun_module
  use flu_binding
  use flu_kinds_module, only: double_k
  use aot_quadruple_top_module, only: quad_k
  use aot_fun_declaration_module, only: aot_fun_type
  use aot_table_module, only: aot_table_from_1Darray

  implicit none

  private

  public :: aot_fun_put

  !> Put an argument into the lua function.
  !!
  !! Arguments have to be in order, first put the first argument then the second
  !! and so on.
  !! Here we add support for quadruple precision numbers
  interface aot_fun_put
    module procedure aot_fun_put_quadruple
    module procedure aot_fun_put_quadruple_v
  end interface aot_fun_put

contains

  !> Put an argument of type double into the list of arguments for the function.
  subroutine aot_fun_put_quadruple(L, fun, arg)
    type(flu_state) :: L !< Handle for the Lua script.

    !> Handle of the function, this argument should be put into.
    type(aot_fun_type) :: fun

    !> Actual argument to hand over to the Lua function.
    real(kind=quad_k), intent(in) :: arg

    real(kind=double_k) :: locarg

    ! Only do something, if the function is actually properly defined.
    if (fun%handle /= 0) then

      locarg = real(arg, kind=double_k)

      ! If the function was executed before this call, it has to be
      ! reset.
      if (fun%arg_count == -1) then
        ! Set the top of the stack to the reference of the function.
        ! Discarding anything above it.
        call flu_settop(L, fun%handle)
        ! Push a copy of the function itself on the stack again, before
        ! adding arguments, to savely survive popping of the function
        ! upon execution.
        call flu_pushvalue(L, fun%handle)
        ! Increase the argument count to 0 again (really start counting
        ! arguments afterwards.
        fun%arg_count = fun%arg_count+1
      end if

      call flu_pushNumber(L, locarg)
      fun%arg_count = fun%arg_count+1
    end if

  end subroutine aot_fun_put_quadruple

  !> Put an array of quadruples into the list of arguments for the
  !! function.
  subroutine aot_fun_put_quadruple_v(L, fun, arg)
    type(flu_state) :: L !< Handle for the Lua script.

    !> Handle of the function, this argument should be put into.
    type(aot_fun_type) :: fun

    !> Actual argument to hand over to the Lua function.
    real(kind=quad_k), intent(in) :: arg(:)

    real(kind=double_k) :: locarg(size(arg))
    integer :: thandle

    ! Only do something, if the function is actually properly defined.
    if (fun%handle /= 0) then

      locarg = real(arg, kind=double_k)

      ! If the function was executed before this call, it has to be
      ! reset.
      if (fun%arg_count == -1) then
        ! Set the top of the stack to the reference of the function.
        ! Discarding anything above it.
        call flu_settop(L, fun%handle)
        ! Push a copy of the function itself on the stack again, before
        ! adding arguments, to savely survive popping of the function
        ! upon execution.
        call flu_pushvalue(L, fun%handle)
        ! Increase the argument count to 0 again (really start counting
        ! arguments afterwards.
        fun%arg_count = fun%arg_count+1
      end if

      call aot_table_from_1Darray(L, thandle, locarg)
      fun%arg_count = fun%arg_count+1
    end if

  end subroutine aot_fun_put_quadruple_v

end module aot_quadruple_fun_module

max_stars_repo_name: zerothi/aotus

max_issues_repo_name: zerothi/aotus

max_forks_repo_name: zerothi/aotus


text:
      subroutine clawpack5_rpn2cons_manifold(ixy,maxm,meqn,
     &         mwaves,maux,mbc, mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
      implicit none

      integer ixy, maxm, meqn, mwaves,maux, mbc, mx

      double precision wave(meqn, mwaves, 1-mbc:maxm+mbc)
      double precision    s(mwaves, 1-mbc:maxm+mbc)
      double precision   ql(meqn, 1-mbc:maxm+mbc)
      double precision   qr(meqn, 1-mbc:maxm+mbc)
      double precision amdq(meqn, 1-mbc:maxm+mbc)
      double precision apdq(meqn, 1-mbc:maxm+mbc)
      double precision auxl(maux, 1-mbc:maxm+mbc)
      double precision auxr(maux, 1-mbc:maxm+mbc)

      integer :: icom,jcom
      double precision :: dtcom,dxcom,dycom,tcom
      common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom

      integer i, idir
      double precision qll,qrr
      double precision urrot, ulrot, g, uhat

      idir = ixy-1
      do i = 2-mbc, mx+mbc
         g = auxl(6+idir,i)  !! Edge length

c        # left-right : 2,3
c        # top-bottom : 4,5
         urrot = g*auxl(2 + 2*idir,i)   !! Left edge of right cell
         ulrot = g*auxl(3 + 2*idir,i-1)   !! Right edge of left cell

         qrr = ql(1,i)
         qll = qr(1,i-1)

c        # Use Roe-average values
         uhat = (ulrot + urrot)/2.d0

         if (uhat .ge. 0) then
            amdq(1,i) = 0.d0
            apdq(1,i) = urrot*qrr - ulrot*qll
         else
            amdq(1,i) = urrot*qrr - ulrot*qll
            apdq(1,i) = 0.d0
         endif
         wave(1,1,i) = urrot*qrr - ulrot*qll
         s(1,i) = uhat
      enddo

      return
      end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
-makelib ies/xil_defaultlib -sv \
  "C:/Xilinx/Vivado/2016.2/data/ip/xpm/xpm_cdc/hdl/xpm_cdc.sv" \
-endlib
-makelib ies/xpm \
  "C:/Xilinx/Vivado/2016.2/data/ip/xpm/xpm_VCOMP.vhd" \
-endlib
-makelib ies/dist_mem_gen_v8_0_10 \
  "../../../ipstatic/dist_mem_gen_v8_0_10/simulation/dist_mem_gen_v8_0.v" \
-endlib
-makelib ies/xil_defaultlib \
  "../../../../Pipeline.srcs/sources_1/ip/instructionROM/sim/instructionROM.v" \
-endlib
-makelib ies/xil_defaultlib \
  glbl.v
-endlib

max_stars_repo_name: HandsomeBrotherShuaiLi/pipelineCPU

max_issues_repo_name: HandsomeBrotherShuaiLi/pipelineCPU

max_forks_repo_name: HandsomeBrotherShuaiLi/pipelineCPU


text:
      REAL FUNCTION  cozena( day, hour, dlat, dlon )
C+---------------------------------------------------------------------+
C!       Reference:    Woolf, H.M., NASA TM X-1646.                    !
C!                     On the computation of solar elevation           !
C!                     angles and the determination of sunrise         !
C!                     and sunset times.                               !
C+---------------------------------------------------------------------+
C!       D a t a               v a r i a b l e s:                      !
C+---------------------------------------------------------------------+
C!       epsiln   eccentricity of earth orbit                          !
C!       sinob    sine of obliquity of ecliptic                        !
C!       dpy      days per year (365.242)                              !
C!       dph      degrees per hour (360./24)                           !
C+---------------------------------------------------------------------+
C!       I n p u t             v a r i a b l e s:                      !
C+---------------------------------------------------------------------+
C!       day      days of year (jan. 1 is day 1)                       !
C!       hour     hours of greenwich time                              !
C!       dlon     degree longitude of grid point,                      !
C!                counted positive west of greenwich                   !
C!       dlat     degree latitude  of grid point,                      !
C!                counted positive northern hemisphere                 !
C+---------------------------------------------------------------------+
C!       O u t p u t           v a r i a b l e s:                      !
C+---------------------------------------------------------------------+
C!       cozena   cosine of zenith angle   ( reference -- eq. 1.1 )    !
C+---------------------------------------------------------------------+
C!       I n t e r n a l       v a r i a b l e s:                      !
C+---------------------------------------------------------------------+
C!       dpr      degree/radian (*57.29578)                            !
C!       rpd      radian/degree (*.01745329)                           !
C!       dang     angle measured from perihelion, in radians           !
C!                which is taken as midnight jan. 1                    !
C!       homp     hours of meridian passage or true solar noon         !
C!                ( reference -- eq. 1.6 )                             !
C!       hang     hour angle, a measure of the longitudinal distance   !
C!                to the sun from the point calculated  ( eq. 1.5 )    !
C!       sindlt   sine of declination angle  ( eq. 1.2 )               !
C!       cosdlt   cosine of declination angle                          !
C!       sigma    reference -- eq. 1.3a                                !
C!       ang      reference -- eq. 1.3b                                !
C+---------------------------------------------------------------------+
      REAL ang, cosdlt, dang, day, dlat,
     $     dlon, dph, dpy, epsiln,
     $     hang, homp, hour, pi, rpd, dpr, sigma, sindlt, sinob
      DATA   epsiln/.016733/, sinob/.3978/, dpy/365.242/,
     $       dph/15.0/, pi/3.1415926535898/
      rpd    = pi/180.0
      dpr    = 1.0/rpd
      dang   = 2.0*pi*(day - 1.0)/dpy
      homp   = 12.0 + 0.123570*SIN( dang )
     $     - 0.004289*COS( dang )
     $     + 0.153809*SIN( 2.0*dang )
     $     + 0.060783*COS( 2.0*dang )
      hang   = dph*( hour-homp ) - dlon
      ang    = 279.9348*rpd + dang
      sigma  = ( ang*dpr + 0.4087*SIN( ang )
     $     + 1.8724*COS( ang )
     $     - 0.0182*SIN( 2.0*ang )
     $     + 0.0083*COS( 2.0*ang ) )*rpd
      sindlt = sinob*SIN( sigma )
      cosdlt = SQRT( 1.0-sindlt**2 )
      cozena = sindlt*SIN( rpd*dlat ) +
     $         cosdlt*COS( rpd*dlat )*COS( rpd*hang )
      RETURN
      END


!*==zpttrs.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b ZPTTRS
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download ZPTTRS + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE ZPTTRS( UPLO, N, NRHS, D, E, B, LDB, INFO )
!
!       .. Scalar Arguments ..
!       CHARACTER          UPLO
!       INTEGER            INFO, LDB, N, NRHS
!       ..
!       .. Array Arguments ..
!       DOUBLE PRECISION   D( * )
!       COMPLEX*16         B( LDB, * ), E( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> ZPTTRS solves a tridiagonal system of the form
!>    A * X = B
!> using the factorization A = U**H *D* U or A = L*D*L**H computed by ZPTTRF.
!> D is a diagonal matrix specified in the vector D, U (or L) is a unit
!> bidiagonal matrix whose superdiagonal (subdiagonal) is specified in
!> the vector E, and X and B are N by NRHS matrices.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] UPLO
!> \verbatim
!>          UPLO is CHARACTER*1
!>          Specifies the form of the factorization and whether the
!>          vector E is the superdiagonal of the upper bidiagonal factor
!>          U or the subdiagonal of the lower bidiagonal factor L.
!>          = 'U':  A = U**H *D*U, E is the superdiagonal of U
!>          = 'L':  A = L*D*L**H, E is the subdiagonal of L
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the tridiagonal matrix A.  N >= 0.
!> \endverbatim
!>
!> \param[in] NRHS
!> \verbatim
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of columns
!>          of the matrix B.  NRHS >= 0.
!> \endverbatim
!>
!> \param[in] D
!> \verbatim
!>          D is DOUBLE PRECISION array, dimension (N)
!>          The n diagonal elements of the diagonal matrix D from the
!>          factorization A = U**H *D*U or A = L*D*L**H.
!> \endverbatim
!>
!> \param[in] E
!> \verbatim
!>          E is COMPLEX*16 array, dimension (N-1)
!>          If UPLO = 'U', the (n-1) superdiagonal elements of the unit
!>          bidiagonal factor U from the factorization A = U**H*D*U.
!>          If UPLO = 'L', the (n-1) subdiagonal elements of the unit
!>          bidiagonal factor L from the factorization A = L*D*L**H.
!> \endverbatim
!>
!> \param[in,out] B
!> \verbatim
!>          B is COMPLEX*16 array, dimension (LDB,NRHS)
!>          On entry, the right hand side vectors B for the system of
!>          linear equations.
!>          On exit, the solution vectors, X.
!> \endverbatim
!>
!> \param[in] LDB
!> \verbatim
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: if INFO = -k, the k-th argument had an illegal value
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date June 2016
!
!> \ingroup complex16PTcomputational
!
!  =====================================================================
      SUBROUTINE ZPTTRS(Uplo,N,Nrhs,D,E,B,Ldb,Info)
      IMPLICIT NONE
!*--ZPTTRS125
!
!  -- LAPACK computational routine (version 3.7.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     June 2016
!
!     .. Scalar Arguments ..
      CHARACTER Uplo
      INTEGER Info , Ldb , N , Nrhs
!     ..
!     .. Array Arguments ..
      DOUBLE PRECISION D(*)
      COMPLEX*16 B(Ldb,*) , E(*)
!     ..
!
!  =====================================================================
!
!     .. Local Scalars ..
      LOGICAL upper
      INTEGER iuplo , j , jb , nb
!     ..
!     .. External Functions ..
      INTEGER ILAENV
      EXTERNAL ILAENV
!     ..
!     .. External Subroutines ..
      EXTERNAL XERBLA , ZPTTS2
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC MAX , MIN
!     ..
!     .. Executable Statements ..
!
!     Test the input arguments.
!
      Info = 0
      upper = (Uplo=='U' .OR. Uplo=='u')
      IF ( .NOT.upper .AND. .NOT.(Uplo=='L' .OR. Uplo=='l') ) THEN
         Info = -1
      ELSEIF ( N<0 ) THEN
         Info = -2
      ELSEIF ( Nrhs<0 ) THEN
         Info = -3
      ELSEIF ( Ldb<MAX(1,N) ) THEN
         Info = -7
      ENDIF
      IF ( Info/=0 ) THEN
         CALL XERBLA('ZPTTRS',-Info)
         RETURN
      ENDIF
!
!     Quick return if possible
!
      IF ( N==0 .OR. Nrhs==0 ) RETURN
!
!     Determine the number of right-hand sides to solve at a time.
!
      IF ( Nrhs==1 ) THEN
         nb = 1
      ELSE
         nb = MAX(1,ILAENV(1,'ZPTTRS',Uplo,N,Nrhs,-1,-1))
      ENDIF
!
!     Decode UPLO
!
      IF ( upper ) THEN
         iuplo = 1
      ELSE
         iuplo = 0
      ENDIF
!
      IF ( nb>=Nrhs ) THEN
         CALL ZPTTS2(iuplo,N,Nrhs,D,E,B,Ldb)
      ELSE
         DO j = 1 , Nrhs , nb
            jb = MIN(Nrhs-j+1,nb)
            CALL ZPTTS2(iuplo,N,jb,D,E,B(1,j),Ldb)
         ENDDO
      ENDIF
!
!
!     End of ZPTTRS
!
      END SUBROUTINE ZPTTRS


text:
!##############################################################################
!# ****************************************************************************
!# <name> interpolation </name>
!# ****************************************************************************
!#
!# <purpose>
!#   This module contains some interpolation routines
!# </purpose>
!##############################################################################

module interpolation

  use fsystem
  use triangulation

implicit none

contains

!<subroutine>

  subroutine Interpol_Press(DpEl,DpMid,rtriangulation,ipar)

  !<description>
    ! This subroutine enables the interpolation of the pressure
    ! from piecewise constant approximation to piecewice linear
    ! and back, depending on the ipar value
    ! ipar == 0 : convert piecewise constant DPC to piecewise linear DPL
    ! ipar == 1 : convert piecewise linear DPL to piecewise constant DPC
  !</description>

  !<input>
    type (t_triangulation), intent(IN) :: rtriangulation
    real(DP), intent(INOUT), dimension(:) :: DpEl,DpMid
    integer(I32),intent(IN) ::  ipar
  !</input>

  !<output>

  !</output>

!</subroutine>

  ! local variables
  integer(I32), dimension(:,:), pointer :: p_Kmid, p_Kadj

  integer(I32) :: iel,iadj,imid,ive
  real(DP) :: dph

  call storage_getbase_int2D(rtriangulation%h_IedgesAtElement,p_Kmid)
  call storage_getbase_int2D(rtriangulation%h_IneighboursAtElement,p_Kadj)

  if (ipar == 0) then

    do iel=1,rtriangulation%NEL
      dph = DpEl(iel)

      do ive=1,4
        iadj=p_Kadj(ive,iel)
        imid=p_Kmid(ive,iel)-rtriangulation%NVT

        if (iadj == 0) DpMid(imid)=dph
        if (iadj >iel) DpMid(imid)=0.5d0*(dph+DpEl(iadj))

      end do
    end do

  else

    do iel=1,rtriangulation%NEL
      DpEl(iel)=0.25d0*(DpMid(p_Kmid(1,iel)-rtriangulation%NVT)+ &
                        DpMid(p_Kmid(2,iel)-rtriangulation%NVT)+ &
                        DpMid(p_Kmid(3,iel)-rtriangulation%NVT)+ &
                        DpMid(p_Kmid(4,iel)-rtriangulation%NVT))
    end do

  endif

  end subroutine Interpol_Press

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
Pants are Optional is Alex Olenicks poor excuse for a comic published in The California Aggie. Olenick probably got the spot because he is a Alpha Epsilon Pi fraternity brother of Daniel Stone Daniel Stones. Users/SamSellers Sam Sellers, the guy who draws the Aggie comic Paperclips, http://www.livejournal.com/users/time_master/32305.html took a dump on this comic on Nov. 30 2005, but no one really got it.

Pants are Optional used to be found http://www.californiaaggie.com/features/comics/?id222 here.

 Why this comic is good
  Well..... its not all THAT bad. Users/PaulThober PT

 Why this comic sucks
  The artistic style is stolen from a popular webcomic, http://pennyarcade.com Penny Arcade
    Hear hear!  Users/JesseSingh js
        I know someone who complained about that to Daniel Stone. He said hed talk to him, but Ill have to ask him if he got a reply from Stone yet.  11/4 Comic in the Aggie vs http://www.pennyarcade.com/comic/2000/10/18/ this one.  11/4 comic in the Aggie is not online....?
  The name of the strip is stole from a popular webcomic, http://ctrlaltdelonline.com/ Ctrl+Alt+Del, where pants are optional
  The comic often has no discernible punchline
  The art is not only a knockoff, but is poor and ignores subtle details such as texture, shading and backgrounds.
  The strip is scanned at 72 DPI and/or saved in a lossy format. Not only does the art suffer, but the text is often difficult, and occasionally impossible, to read.
  Rather than letter by hand, Olenick uses wiki:WikiPedia:Comic_Sans Comic Sans, quite possibly the most hideous font ever created by man.

About a very, very similar comic strip idea:  11/4 Comic in the Aggie vs http://www.pennyarcade.com/comic/2000/10/18/ this one.  For some reason, that specific comic was not in the Aggie link above (currently broken), though the one dated the week before was and the one the week after was. Strange, eh? Must have taken it out.

A message was sent 114 to Daniel Stone about the two comics. No reply.   On 119, he sent the message to 6 or 7 editors emails. He got a reply immediately from Stone, about two hours after.

Thank you for your comments about The California Aggie. I did receive your concern on Nov. 4 and I do appologize for not responding sooner. Regarding the comic strip, I have asked the comic artist to come in to discuss the strip. I have looked at the strips you provided and do see a remarkable similarity, though I will give him the opportunity to defend himself.

Again, thank you for your concerns

Daniel Stone Editor in Chief The California Aggie

Then..nothing. A reply was sent on Nov 29th asking about any sort of followup since then.  As of December 3rd, there was no response.  Ive yet to ask again, but I personally doubt he got one in the last three days.

I asked again.  He sent Daniel Stone yet another email asking for a followup December 10th.  Reply came Dec 12, a month after the first one.  Thank you for your comments regarding Alex Olenicks comic in The Aggie. All matters of personnel are kept confidential and are not discussed outside of the management staff. I appreciate your concerns and willingness to bring pertinent information to my attention Daniel Stone

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      LOGICAL FUNCTION CSSIDS (IDSTR)

C  This routine checks the actual set of structure identifiers in
C  the CSS against an expected list, encoded in idstr.  If equal it
C  returns true, else false.

      INTEGER    ARRSIZ, IDUM, IX, ERRIND
      PARAMETER (ARRSIZ = 300)

      INTEGER    EXPIDS(ARRSIZ), EXPSIZ, ACTIDS(ARRSIZ), ACTSIZ
      CHARACTER  IDSTR*(*)
      LOGICAL    SETEQ

      CALL SETVS (IDSTR, EXPIDS, EXPSIZ)
      IF (EXPSIZ .GT. ARRSIZ) THEN
         CALL UNMSG ('Expected size exceeds limits of CSSIDS.')
      ENDIF

C  <inquire structure identifiers> to determine actual list of structures
      CALL PQSID (0, ERRIND, ACTSIZ, IDUM)
      CALL CHKINQ ('pqsid', ERRIND)
      IF (ACTSIZ .NE. EXPSIZ) THEN
         CSSIDS = .FALSE.
         RETURN
      ENDIF
      IF (ACTSIZ .GT. ARRSIZ) THEN
         CALL UNMSG ('Actual size exceeds limits of CSSIDS.')
      ENDIF

      DO 100 IX = 1,ACTSIZ
         CALL PQSID (IX, ERRIND, IDUM, ACTIDS(IX))
         CALL CHKINQ ('pqsid', ERRIND)
100   CONTINUE

      CSSIDS = SETEQ (ACTSIZ, ACTIDS, EXPIDS)

      END


text:
      program w13

C     Derived from w12. ELSEIF added to obtain unfeasible paths

C     Assignments in body: check that the initial condition I=0 and the
C     body exit condition 0<=I<=1 are merged as entry condition for the
C     loop body before the fixpoint is computed.

C     To obtain a precise loop transformer, the transformer for the
C     entered loop and the transformer for the non-entered loop should
C     be computed and their convex hull taken, but it is not useful,
C     since the operation is performed when preconditions are computed.

      integer i

      i = 0

      do while(x.gt.0.)
         if(i.eq.1) then
            i = 0
         elseif(i.eq.0)then
            i = 1
         endif
         read *, x
      enddo

      print *, i

      end


text:
MODULE m_eop_igu

! ----------------------------------------------------------------------
! MODULE: m_eop_igu.f03
! ----------------------------------------------------------------------
! Purpose:
!  Module for calling the modified eop_igu subroutine
!
! ----------------------------------------------------------------------
! Author :	Dr. Thomas Papanikolaou, Geoscience Australia
! Created:	28 August 2018
! ----------------------------------------------------------------------

      IMPLICIT NONE
      !SAVE

Contains

SUBROUTINE eop_igu (mjd, ERP_fname, ERP_days, EOP_days, EOP_Nint, EOP_int)

! ----------------------------------------------------------------------
! Subroutine:  eop_igu.f90
! ----------------------------------------------------------------------
! Purpose:
!  EOP data reading and processing by using:
!  - ERP (Earth Rotation Parameters) data from the ultra-rapid products
!    provided by the IGS (International GNSS Service).
!  - Corrections to the precession-nutation model are obtained from the
!    daily solutions (finals2000A.daily) provided by the International
!    Earth Rotation Service and Reference Systems (IERS)
!    Rapid Service/Prediction Center (RS/PC)
! ----------------------------------------------------------------------
! Input arguments:
! - mjd:			Modified Julian Day number at the required epoch
!					(including fraction of the day)
! - ERP_fname:		IGS ultra-rapid ERP data file name e.g. igu18861_00.erp
! - EOP_days:		EOP data array of the days (data points aplied for interpolation) based on IERS RS/PC EOP data
! - EOP_Nint:           Number of points to use for interpolation
! - arc_length:         Number of hours of arc to read from file
!
! Output arguments:
! - eop_int:		EOP data array at the input epoch
!   				eop_int = [MJD xp yp UT1_UTC LOD dX dY]
!   				MJD:     MJD at the input epoch (including fraction of the day)
!   				x,y:     Polar motion coordinates (arcsec)
!   				UT1_UTC: Difference between UT1 and UTC (sec)
!					dX,dY:   Corrections to Precession-Nutation model (arcsec)
! ----------------------------------------------------------------------
! Dr. Thomas Papanikolaou, Geoscience Australia               March 2016
! ----------------------------------------------------------------------

      USE mdl_precision
      USE mdl_num
      use pod_yaml
      IMPLICIT NONE

! ----------------------------------------------------------------------
! Dummy arguments declaration
! ----------------------------------------------------------------------
! IN
      REAL (KIND = prec_d), INTENT(IN) :: mjd
      CHARACTER (LEN=512), INTENT(IN) :: ERP_fname
      REAL (KIND = prec_d), INTENT(IN), DIMENSION(:,:), ALLOCATABLE :: ERP_days, EOP_days
      INTEGER, INTENT (IN) :: EOP_Nint
! OUT
      REAL (KIND = prec_d), INTENT(OUT) :: EOP_int(EOP_MAX_ARRAY)
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! Local variables declaration
! ----------------------------------------------------------------------
      REAL (KIND = prec_d) :: ERP_igu_data(EOP_Nint,EOP_MAX_ARRAY), ERP_int(EOP_MAX_ARRAY), mjd_t
      LOGICAL :: igu_flag
      INTEGER (KIND = prec_int8) :: mjd_UTC_day
      REAL (KIND = prec_d) :: mjd_ar(EOP_Nint), Xpole_ar(EOP_Nint), Ypole_ar(EOP_Nint), UT1UTC_ar(EOP_Nint), LOD_ar(EOP_Nint)
      REAL (KIND = prec_d) :: Xerr_ar(EOP_Nint), Yerr_ar(EOP_Nint), UT1err_ar(EOP_Nint), LODerr_ar(EOP_Nint)
      REAL (KIND = prec_d) :: mjd_int, Xpole_int, Ypole_int, UT1UTC_int, LOD_int
      REAL (KIND = prec_d) :: Xerr_int, Yerr_int, UT1err_int, LODerr_int
      REAL (KIND = prec_d) :: dX_eop, dY_eop
      INTEGER (KIND = prec_int4) :: i, j, sz1_EOP, sz2_EOP, EOP_nint_used, saved_row_count = 0
      INTEGER (KIND = prec_int4) :: lo, hi, offset
      INTEGER MAX_EOP_DATA
      parameter (MAX_EOP_DATA = 30)
      REAL (KIND = prec_d) :: ERP_read_data(MAX_EOP_DATA, EOP_MAX_ARRAY)

      igu_flag = .false.
      ERP_igu_data = 0.d0

! ----------------------------------------------------------------------
! check we don't have the required data already
      if (allocated(ERP_days)) then
          saved_row_count = SIZE(ERP_days, DIM=1)
          do i = 1, saved_row_count
              if (ERP_days(i,EOP_MJD) == 0.d0) then
                   saved_row_count = i-1
                   exit
              end if
              ERP_read_data(i,:) = ERP_days(i, :)
          end do
         DO i = 1, saved_row_count - 1
         if ((ERP_read_data(i, EOP_MJD) <= mjd) .and. (mjd <= ERP_read_data(i+1, EOP_MJD))) then
            igu_flag = .true.
            if (i > (EOP_Nint / 2) .and. ((i + 1) <= (saved_row_count - EOP_Nint / 2))) then
               if (MOD(EOP_Nint, 2) == 1) then
                   lo = i-EOP_Nint/2
                   hi = i+1+EOP_Nint/2
                   if (lo == 0) then
                       lo = lo + 1
                       hi = hi + 1
                   end if
               else
                   lo = i-EOP_Nint/2+1
                   hi = i+EOP_Nint/2
               end if
               ERP_igu_data = ERP_read_data(lo:hi, :)
               EOP_Nint_used = EOP_Nint
            else if (i <= EOP_nint/2) then
                ERP_igu_data = ERP_read_data(1:EOP_Nint,:)
                EOP_Nint_used = EOP_nint
            else if (i >= (saved_row_count - EOP_Nint/2)) then
                ERP_igu_data = ERP_read_data(saved_row_count - EOP_Nint + 1:saved_row_count, :)
                EOP_nint_used = EOP_Nint
            end if
            exit
         end if
         end do
         if (.not. igu_flag) then
             !print *, "mjd=", mjd, ", last row of data is mjd ", ERP_read_data(saved_row_count, EOP_MJD)
             if (mjd < ERP_read_data(saved_row_count, EOP_MJD) + 1.0d0) then
                 igu_flag = .true.
                 ERP_igu_data = ERP_read_data(saved_row_count - EOP_Nint + 1:saved_row_count, :)
                 EOP_Nint_used = EOP_Nint
             else if (mjd < ERP_read_data(1, EOP_MJD)) then
                 igu_flag = .true.
                 ERP_igu_data = ERP_read_data(1:EOP_nint,:)
                 EOP_Nint_used = EOP_Nint
             end if
         end if
     end if

! ----------------------------------------------------------------------
! ERP data reading
      if (.not. igu_flag) then
          if (MOD(EOP_Nint, 2) == 1) then
              !EOP_Nint is odd, use Nint/2 points below if mjd is above 0.5, else Nint/2 + 1 points below
              mjd_t = mjd - int(mjd)
              if (mjd_t > 0.5) then
                  lo = int(mjd) - EOP_Nint/2 + 1
                  hi = INT(mjd) + EOP_Nint/2 + 1
              else
                  lo = int(mjd) - EOP_Nint/2
                  hi = int(mjd) + EOP_Nint/2
              end if
          else
              !EOP_nint is even use half above and half below precisely
              lo = int(mjd) - EOP_Nint/2 + 1
              hi = int(mjd) + EOP_Nint/2
          end if
          offset = 1
          EOP_nint_used = 0
          do i = hi, lo, -1
              mjd_t = float(i)
              CALL erp_igu (ERP_fname, mjd_t, ERP_int, igu_flag)
              if (igu_flag) then
                 erp_igu_data (i - lo + offset, :) = ERP_int
                 if (EOP_Nint_used == 0) then
                    EOP_Nint_used = i - lo + offset
                 end if
              else
                offset = offset + 1
              end if
          end do
      end if
      if (.not. igu_flag .and. offset > 2 ) then
         PRINT *,"--------------------------------------------------------"
         PRINT *, "Warning error: Subroutine erp_igu.f90"
         PRINT *, "Input epoch is out of the range covered by the IGS ultra-rapid ERP file"
         PRINT *, "Check the input ", ERP_fname
         PRINT *,"--------------------------------------------------------"
         STOP  ! END PROGRAM
      end if
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! ERP interpolation
      ERP_int = 0.d0

      mjd_int = mjd
      mjd_ar    = ERP_igu_data(1:EOP_Nint_used,EOP_MJD)
      Xpole_ar  = ERP_igu_data(1:EOP_Nint_used,EOP_X)
      Ypole_ar  = ERP_igu_data(1:EOP_Nint_used,EOP_Y)
      UT1UTC_ar = ERP_igu_data(1:EOP_Nint_used,EOP_UT1)
      LOD_ar    = ERP_igu_data(1:EOP_Nint_used,EOP_LOD)
      Xerr_ar = ERP_igu_data(1:EOP_Nint_used,EOP_X_ERR)
      Yerr_ar = ERP_igu_data(1:EOP_Nint_used,EOP_Y_ERR)
      UT1err_ar = ERP_igu_data(1:EOP_Nint_used,EOP_UT1_ERR)
      LODerr_ar = ERP_igu_data(1:EOP_Nint_used,EOP_LOD_ERR)
      ! do we need to correct, or should that have already been applied?
      call interp_iers(mjd_ar, Xpole_ar, Ypole_ar, UT1UTC_ar, EOP_Nint_used, EOP_Nint_used, &
              mjd_int, Xpole_int, Ypole_int, UT1UTC_int)
      !CALL LAGINT(mjd_ar, Xpole_ar, EOP_Nint_used, EOP_Nint_used, mjd_int, Xpole_int)
      !CALL LAGINT(mjd_ar, Ypole_ar, EOP_Nint_used, EOP_Nint_used, mjd_int, Ypole_int)
      !CALL LAGINT(mjd_ar, UT1UTC_ar, EOP_Nint_used, EOP_Nint_used, mjd_int, UT1UTC_int)
      CALL LAGINT(mjd_ar, LOD_ar, EOP_Nint_used, EOP_Nint_used, mjd_int, LOD_int)
      CALL LAGINT(mjd_ar, Xerr_ar, EOP_Nint_used, EOP_Nint_used, mjd_int, Xerr_int)
      CALL LAGINT(mjd_ar, Yerr_ar, EOP_Nint_used, EOP_Nint_used, mjd_int, Yerr_int)
      CALL LAGINT(mjd_ar, UT1err_ar, EOP_Nint_used, EOP_Nint_used, mjd_int, UT1err_int)
      CALL LAGINT(mjd_ar, LODerr_ar, EOP_Nint_used, EOP_Nint_used, mjd_int, LODerr_int)

      ERP_int (EOP_MJD) = mjd_int
      ERP_int (EOP_X) = Xpole_int
      ERP_int (EOP_Y) = Ypole_int
      ERP_int (EOP_UT1) = UT1UTC_int
      ERP_int (EOP_LOD) = LOD_int
      ERP_int (EOP_X_ERR) = Xerr_int
      ERP_int (EOP_Y_ERR) = Yerr_int
      ERP_int (EOP_UT1_ERR) = ut1err_int
      ERP_int (EOP_LOD_ERR) = LODerr_int
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
! dX,dY : Corrections w.r.t Precession-Nutation model
      !mjd_UTC_day = INT (mjd)
      !CALL eop_finals2000A (EOP_fname, mjd_UTC_day , EOP_data)
      !dX = EOP_data(6)
      !dY = EOP_data(7)
      !----------------------------------------------------------
      ! init next two var to something sensible
dX_eop = 0.d0
dY_eop = 0.d0
sz1_EOP = SIZE (EOP_days,DIM=1)
sz2_EOP = SIZE (EOP_days,DIM=2)
DO i = 1 , sz1_EOP
If (mjd_int == EOP_days(i,EOP_MJD) ) then
! dX,dY (arcsec)
      dX_eop = EOP_days(i,EOP_DX)
      dY_eop = EOP_days(i,EOP_DY)
End If
END DO
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
      EOP_int (EOP_MJD) = mjd_int
      EOP_int (EOP_X) = Xpole_int
      EOP_int (EOP_Y) = Ypole_int
      EOP_int (EOP_UT1) = UT1UTC_int
      EOP_int (EOP_LOD) = LOD_int
	  ! dX,dY : Precession-Nutation model corrections
      EOP_int (EOP_DX) = dX_eop !EOP_data (6)
      EOP_int (EOP_DY) = dY_eop !EOP_data (7)
! ----------------------------------------------------------------------
      EOP_int (EOP_X_ERR) = Xerr_int
      EOP_int (EOP_Y_ERR) = Yerr_int
      EOP_int (EOP_UT1_ERR) = ut1err_int
      EOP_int (EOP_LOD_ERR) = LODerr_int

!      PRINT *,"--------------------------------------------------------"
!      PRINT *, "ERP_igu_data"
!      PRINT *, ERP_igu_data
!      PRINT *, "ERP_int"
!      PRINT *, ERP_int
!      PRINT *,"--------------------------------------------------------"

END SUBROUTINE

END

max_stars_repo_name: HiTMonitor/ginan

max_issues_repo_name: HiTMonitor/ginan

max_forks_repo_name: HiTMonitor/ginan


text:
!
! Tioga is a library for overset grid assembly on parallel distributed systems
! Copyright (C) 2015 Jay Sitaraman
!
subroutine getSpherePart(myid,numprocs,arange)
!
integer, intent(in) :: myid
integer, intent(in) :: numprocs
real*8, intent(out) :: arange(4)
!
integer :: i,nc,j,k
real*8 :: Ap,thetap,delta,Ai,a,thpolar
real*8, allocatable :: y(:)
integer, allocatable :: m(:)
real*8 :: pi
!
! begin
!
pi=acos(-1.)
Ap=4*pi/numprocs
!
thetap=asin(1-2.0/numprocs)
thpolar=pi*0.5-thetap
!
nc=nint((pi-2*thpolar)/sqrt(Ap))
delta=(pi-2*thpolar)/nc
!
allocate(y(nc),m(nc))
!
do i=1,nc
   Ai=2*pi*(sin(thetap-(i-1)*delta)-sin(thetap-i*delta))
   y(i)=Ai/Ap
enddo
!
a=0
do i=1,nc
   m(i)=nint(y(i)+a)
   a=a+(y(i)-m(i))
enddo
!
k=1
if (myid==k) then
  arange(1)=thetap
  arange(2)=pi*0.5
  arange(3)=0.
  arange(4)=2*pi
endif
!
do i=1,nc
   do j=1,m(i)
     k=k+1
     if (myid==k) then
        arange(1)=thetap-i*delta
	arange(2)=arange(1)+delta
        arange(3)=(j-1)*2*pi/m(i)
        arange(4)=arange(3)+2*pi/m(i)
	return
     endif
   enddo
enddo
!
k=numprocs
if (myid==k) then
  arange(1)=-pi*0.5
  arange(2)=-thetap
  arange(3)=0.
  arange(4)=2*pi
endif
!
!
return
end subroutine getSpherePart

max_stars_repo_name: jsitaraman/fvsand

max_issues_repo_name: jsitaraman/fvsand

max_forks_repo_name: jsitaraman/fvsand


! -----------------------------------------------------------------------------
! This file was automatically created by SARAH version SARAHVERSION
! SARAH References: arXiv:0806.0538, 0909.2863, 1002.0840, 1207.0906, 1309.7223
! (c) Florian Staub, 2013
! ------------------------------------------------------------------------------
! File created at 13:23 on 7.1.2017
! ----------------------------------------------------------------------

Module RGEs_SM_HC

Use Control
Use Settings
Use Mathematics

Logical,Private,Save::OnlyDiagonal

Real(dp),Parameter::id3R(3,3)=&
   & Reshape(Source=(/&
   & 1,0,0,&
 &0,1,0,&
 &0,0,1&
 &/),shape=(/3,3/))
Contains

Subroutine GToParameters62_SM(g,g1,g2,g3,Lam,Yu,Yd,Ye,Mu,v)

Implicit None
Real(dp), Intent(in) :: g(62)
Real(dp),Intent(out) :: g1,g2,g3,v

Complex(dp),Intent(out) :: Lam,Yu(3,3),Yd(3,3),Ye(3,3),Mu

Integer i1, i2, i3, i4, SumI

Iname = Iname +1
NameOfUnit(Iname) = 'GToParameters62'

g1= g(1)
g2= g(2)
g3= g(3)
Lam= Cmplx(g(4),g(5),dp)
Do i1 = 1,3
Do i2 = 1,3
SumI = (i2-1) + (i1-1)*3
SumI = SumI*2
Yu(i1,i2) = Cmplx( g(SumI+6), g(SumI+7), dp)
End Do
 End Do

Do i1 = 1,3
Do i2 = 1,3
SumI = (i2-1) + (i1-1)*3
SumI = SumI*2
Yd(i1,i2) = Cmplx( g(SumI+24), g(SumI+25), dp)
End Do
 End Do

Do i1 = 1,3
Do i2 = 1,3
SumI = (i2-1) + (i1-1)*3
SumI = SumI*2
Ye(i1,i2) = Cmplx( g(SumI+42), g(SumI+43), dp)
End Do
 End Do

Mu= Cmplx(g(60),g(61),dp)
v= g(62)
Do i1=1,62
If (g(i1).ne.g(i1)) Then
 Write(*,*) "NaN appearing in ",NameOfUnit(Iname)
 Write(*,*) "At position ", i1
 Call TerminateProgram
End if
End do
Iname = Iname - 1

End Subroutine GToParameters62_SM

Subroutine ParametersToG62_SM(g1,g2,g3,Lam,Yu,Yd,Ye,Mu,v,g)

Implicit None
Real(dp), Intent(out) :: g(62)
Real(dp), Intent(in) :: g1,g2,g3,v

Complex(dp), Intent(in) :: Lam,Yu(3,3),Yd(3,3),Ye(3,3),Mu

Integer i1, i2, i3, i4, SumI

Iname = Iname +1
NameOfUnit(Iname) = 'ParametersToG62'

g(1) = g1
g(2) = g2
g(3) = g3
g(4) = Real(Lam,dp)
g(5) = Aimag(Lam)
Do i1 = 1,3
Do i2 = 1,3
SumI = (i2-1) + (i1-1)*3
SumI = SumI*2
g(SumI+6) = Real(Yu(i1,i2), dp)
g(SumI+7) = Aimag(Yu(i1,i2))
End Do
End Do

Do i1 = 1,3
Do i2 = 1,3
SumI = (i2-1) + (i1-1)*3
SumI = SumI*2
g(SumI+24) = Real(Yd(i1,i2), dp)
g(SumI+25) = Aimag(Yd(i1,i2))
End Do
End Do

Do i1 = 1,3
Do i2 = 1,3
SumI = (i2-1) + (i1-1)*3
SumI = SumI*2
g(SumI+42) = Real(Ye(i1,i2), dp)
g(SumI+43) = Aimag(Ye(i1,i2))
End Do
End Do

g(60) = Real(Mu,dp)
g(61) = Aimag(Mu)
g(62) = v
Iname = Iname - 1

End Subroutine ParametersToG62_SM

Subroutine rge62_SM(len, T, GY, F)
Implicit None
Integer, Intent(in) :: len
Real(dp), Intent(in) :: T, GY(len)
Real(dp), Intent(out) :: F(len)
Integer :: i1,i2,i3,i4
Integer :: j1,j2,j3,j4,j5,j6,j7
Real(dp) :: q
Real(dp) :: g1,betag11,betag12,Dg1,g2,betag21,betag22,Dg2,g3,betag31,betag32,         &
& Dg3,v,betav1,betav2,Dv, betag33, betaLam3, BetaYu3, Yt, Z3
Complex(dp) :: Lam,betaLam1,betaLam2,DLam,Yu(3,3),betaYu1(3,3),betaYu2(3,3)           &
& ,DYu(3,3),adjYu(3,3),Yd(3,3),betaYd1(3,3),betaYd2(3,3),DYd(3,3),adjYd(3,3)             &
& ,Ye(3,3),betaYe1(3,3),betaYe2(3,3),DYe(3,3),adjYe(3,3),Mu,betaMu1,betaMu2,DMu
Complex(dp) :: YdadjYd(3,3),YeadjYe(3,3),YuadjYu(3,3),adjYdYd(3,3),adjYeYe(3,3),adjYuYu(3,3),        &
& YdadjYdYd(3,3),YdadjYuYu(3,3),YeadjYeYe(3,3),YuadjYdYd(3,3),YuadjYuYu(3,3),            &
& adjYdYdadjYd(3,3),adjYeYeadjYe(3,3),adjYuYuadjYu(3,3),YdadjYdYdadjYd(3,3),             &
& YeadjYeYeadjYe(3,3),YuadjYuYuadjYu(3,3)

Complex(dp) :: YuadjYd(3,3),adjYuYuadjYd(3,3),YdadjYuYuadjYd(3,3),YuadjYdYdadjYd(3,3),               &
& YuadjYuYuadjYd(3,3),adjYdYdadjYdYd(3,3),adjYdYdadjYuYu(3,3),adjYeYeadjYeYe(3,3),       &
& adjYuYuadjYdYd(3,3),adjYuYuadjYuYu(3,3),YdadjYdYdadjYdYd(3,3),YdadjYdYdadjYuYu(3,3),   &
& YdadjYuYuadjYdYd(3,3),YdadjYuYuadjYuYu(3,3),YeadjYeYeadjYeYe(3,3),YuadjYdYdadjYdYd(3,3),&
& YuadjYdYdadjYuYu(3,3),YuadjYuYuadjYdYd(3,3),YuadjYuYuadjYuYu(3,3),adjYdYdadjYdYdadjYd(3,3),&
& adjYdYdadjYuYuadjYd(3,3),adjYeYeadjYeYeadjYe(3,3),adjYuYuadjYdYdadjYd(3,3),            &
& adjYuYuadjYuYuadjYd(3,3),adjYuYuadjYuYuadjYu(3,3),YdadjYdYdadjYdYdadjYd(3,3),          &
& YdadjYdYdadjYuYuadjYd(3,3),YdadjYuYuadjYdYdadjYd(3,3),YdadjYuYuadjYuYuadjYd(3,3),      &
& YeadjYeYeadjYeYeadjYe(3,3),YuadjYuYuadjYuYuadjYu(3,3)

Complex(dp) :: TrYdadjYd,TrYeadjYe,TrYuadjYu,TrYdadjYdYdadjYd,TrYeadjYeYeadjYe,TrYuadjYuYuadjYu

Complex(dp) :: TrYdadjYuYuadjYd,TrYdadjYdYdadjYdYdadjYd,TrYdadjYdYdadjYuYuadjYd,TrYdadjYuYuadjYdYdadjYd,&
& TrYdadjYuYuadjYuYuadjYd,TrYeadjYeYeadjYeYeadjYe,TrYuadjYuYuadjYuYuadjYu

Real(dp) :: g1p2,g1p3,g1p4,g2p2,g2p3,g2p4,g3p2,g3p3

Complex(dp) :: Lamp2

Real(dp) :: g1p6,g2p6,g3p4

Complex(dp) :: Xip2,Lamp3

Iname = Iname +1
NameOfUnit(Iname) = 'rge62'

OnlyDiagonal = .Not.GenerationMixing
q = t

Call GToParameters62_SM(gy,g1,g2,g3,Lam,Yu,Yd,Ye,Mu,v)

Yt= Yu(3,3)
Z3=1.2020569031595942366_dp

Call Adjungate(Yu,adjYu)
Call Adjungate(Yd,adjYd)
Call Adjungate(Ye,adjYe)
 YdadjYd = Matmul(Yd,adjYd)
Forall(i2=1:3)  YdadjYd(i2,i2) =  Real(YdadjYd(i2,i2),dp)
 YeadjYe = Matmul(Ye,adjYe)
Forall(i2=1:3)  YeadjYe(i2,i2) =  Real(YeadjYe(i2,i2),dp)
 YuadjYu = Matmul(Yu,adjYu)
Forall(i2=1:3)  YuadjYu(i2,i2) =  Real(YuadjYu(i2,i2),dp)
 adjYdYd = Matmul(adjYd,Yd)
Forall(i2=1:3)  adjYdYd(i2,i2) =  Real(adjYdYd(i2,i2),dp)
 adjYeYe = Matmul(adjYe,Ye)
Forall(i2=1:3)  adjYeYe(i2,i2) =  Real(adjYeYe(i2,i2),dp)
 adjYuYu = Matmul(adjYu,Yu)
Forall(i2=1:3)  adjYuYu(i2,i2) =  Real(adjYuYu(i2,i2),dp)
 YdadjYdYd = Matmul(Yd,adjYdYd)
 YdadjYuYu = Matmul(Yd,adjYuYu)
 YeadjYeYe = Matmul(Ye,adjYeYe)
 YuadjYdYd = Matmul(Yu,adjYdYd)
 YuadjYuYu = Matmul(Yu,adjYuYu)
 adjYdYdadjYd = Matmul(adjYd,YdadjYd)
 adjYeYeadjYe = Matmul(adjYe,YeadjYe)
 adjYuYuadjYu = Matmul(adjYu,YuadjYu)
 YdadjYdYdadjYd = Matmul(Yd,adjYdYdadjYd)
Forall(i2=1:3)  YdadjYdYdadjYd(i2,i2) =  Real(YdadjYdYdadjYd(i2,i2),dp)
 YeadjYeYeadjYe = Matmul(Ye,adjYeYeadjYe)
Forall(i2=1:3)  YeadjYeYeadjYe(i2,i2) =  Real(YeadjYeYeadjYe(i2,i2),dp)
 YuadjYuYuadjYu = Matmul(Yu,adjYuYuadjYu)
Forall(i2=1:3)  YuadjYuYuadjYu(i2,i2) =  Real(YuadjYuYuadjYu(i2,i2),dp)
 TrYdadjYd = Real(cTrace(YdadjYd),dp)
 TrYeadjYe = Real(cTrace(YeadjYe),dp)
 TrYuadjYu = Real(cTrace(YuadjYu),dp)
 TrYdadjYdYdadjYd = Real(cTrace(YdadjYdYdadjYd),dp)
 TrYeadjYeYeadjYe = Real(cTrace(YeadjYeYeadjYe),dp)
 TrYuadjYuYuadjYu = Real(cTrace(YuadjYuYuadjYu),dp)
 g1p2 =g1**2
 g1p3 =g1**3
 g1p4 =g1**4
 g2p2 =g2**2
 g2p3 =g2**3
 g2p4 =g2**4
 g3p2 =g3**2
 g3p3 =g3**3
 Lamp2 =Lam**2
 g1p6 =g1**6
 g2p6 =g2**6
 g3p4 =g3**4
 Xip2 =Xi**2
 Lamp3 =Lam**3

!  TwoLoopRGE=.false.

If (TwoLoopRGE) Then
 YuadjYd = Matmul(Yu,adjYd)
 adjYuYuadjYd = Matmul(adjYu,YuadjYd)
 YdadjYuYuadjYd = Matmul(Yd,adjYuYuadjYd)
Forall(i2=1:3)  YdadjYuYuadjYd(i2,i2) =  Real(YdadjYuYuadjYd(i2,i2),dp)
 YuadjYdYdadjYd = Matmul(Yu,adjYdYdadjYd)
 YuadjYuYuadjYd = Matmul(Yu,adjYuYuadjYd)
 adjYdYdadjYdYd = Matmul(adjYd,YdadjYdYd)
Forall(i2=1:3)  adjYdYdadjYdYd(i2,i2) =  Real(adjYdYdadjYdYd(i2,i2),dp)
 adjYdYdadjYuYu = Matmul(adjYd,YdadjYuYu)
 adjYeYeadjYeYe = Matmul(adjYe,YeadjYeYe)
Forall(i2=1:3)  adjYeYeadjYeYe(i2,i2) =  Real(adjYeYeadjYeYe(i2,i2),dp)
 adjYuYuadjYdYd = Matmul(adjYu,YuadjYdYd)
 adjYuYuadjYuYu = Matmul(adjYu,YuadjYuYu)
Forall(i2=1:3)  adjYuYuadjYuYu(i2,i2) =  Real(adjYuYuadjYuYu(i2,i2),dp)
 YdadjYdYdadjYdYd = Matmul(Yd,adjYdYdadjYdYd)
 YdadjYdYdadjYuYu = Matmul(Yd,adjYdYdadjYuYu)
 YdadjYuYuadjYdYd = Matmul(Yd,adjYuYuadjYdYd)
 YdadjYuYuadjYuYu = Matmul(Yd,adjYuYuadjYuYu)
 YeadjYeYeadjYeYe = Matmul(Ye,adjYeYeadjYeYe)
 YuadjYdYdadjYdYd = Matmul(Yu,adjYdYdadjYdYd)
 YuadjYdYdadjYuYu = Matmul(Yu,adjYdYdadjYuYu)
 YuadjYuYuadjYdYd = Matmul(Yu,adjYuYuadjYdYd)
 YuadjYuYuadjYuYu = Matmul(Yu,adjYuYuadjYuYu)
 adjYdYdadjYdYdadjYd = Matmul(adjYd,YdadjYdYdadjYd)
 adjYdYdadjYuYuadjYd = Matmul(adjYd,YdadjYuYuadjYd)
 adjYeYeadjYeYeadjYe = Matmul(adjYe,YeadjYeYeadjYe)
 adjYuYuadjYdYdadjYd = Matmul(adjYu,YuadjYdYdadjYd)
 adjYuYuadjYuYuadjYd = Matmul(adjYu,YuadjYuYuadjYd)
 adjYuYuadjYuYuadjYu = Matmul(adjYu,YuadjYuYuadjYu)
 YdadjYdYdadjYdYdadjYd = Matmul(Yd,adjYdYdadjYdYdadjYd)
Forall(i2=1:3)  YdadjYdYdadjYdYdadjYd(i2,i2) =  Real(YdadjYdYdadjYdYdadjYd(i2,i2),dp)
 YdadjYdYdadjYuYuadjYd = Matmul(Yd,adjYdYdadjYuYuadjYd)
 YdadjYuYuadjYdYdadjYd = Matmul(Yd,adjYuYuadjYdYdadjYd)
 YdadjYuYuadjYuYuadjYd = Matmul(Yd,adjYuYuadjYuYuadjYd)
Forall(i2=1:3)  YdadjYuYuadjYuYuadjYd(i2,i2) =  Real(YdadjYuYuadjYuYuadjYd(i2,i2),dp)
 YeadjYeYeadjYeYeadjYe = Matmul(Ye,adjYeYeadjYeYeadjYe)
Forall(i2=1:3)  YeadjYeYeadjYeYeadjYe(i2,i2) =  Real(YeadjYeYeadjYeYeadjYe(i2,i2),dp)
 YuadjYuYuadjYuYuadjYu = Matmul(Yu,adjYuYuadjYuYuadjYu)
Forall(i2=1:3)  YuadjYuYuadjYuYuadjYu(i2,i2) =  Real(YuadjYuYuadjYuYuadjYu(i2,i2),dp)
 TrYdadjYuYuadjYd = Real(cTrace(YdadjYuYuadjYd),dp)
 TrYdadjYdYdadjYdYdadjYd = Real(cTrace(YdadjYdYdadjYdYdadjYd),dp)
 TrYdadjYdYdadjYuYuadjYd = cTrace(YdadjYdYdadjYuYuadjYd)
 TrYdadjYuYuadjYdYdadjYd = cTrace(YdadjYuYuadjYdYdadjYd)
 TrYdadjYuYuadjYuYuadjYd = Real(cTrace(YdadjYuYuadjYuYuadjYd),dp)
 TrYeadjYeYeadjYeYeadjYe = Real(cTrace(YeadjYeYeadjYeYeadjYe),dp)
 TrYuadjYuYuadjYuYuadjYu = Real(cTrace(YuadjYuYuadjYuYuadjYu),dp)
 g1p6 =g1**6
 g2p6 =g2**6
 g3p4 =g3**4
 Xip2 =Xi**2
 Lamp3 =Lam**3
End If

!--------------------
! g1
!--------------------

betag11  = 41._dp*(g1p3)/10._dp

If (TwoLoopRGE) Then
betag12 = (g1p3*(199._dp*(g1p2) + 135._dp*(g2p2) + 440._dp*(g3p2) - 25._dp*(TrYdadjYd) -        &
&  75._dp*(TrYeadjYe) - 85._dp*(TrYuadjYu)))/50._dp

Dg1 = oo16pi2*( betag11 + oo16pi2 * betag12 )

Else
Dg1 = oo16pi2* betag11
End If

!--------------------
! g2
!--------------------

betag21  = -19._dp*(g2p3)/6._dp

If (TwoLoopRGE) Then
betag22 = (g2p3*(27._dp*(g1p2) + 175._dp*(g2p2) + 360._dp*(g3p2) - 45._dp*(TrYdadjYd) -         &
&  15._dp*(TrYeadjYe) - 45._dp*(TrYuadjYu)))/30._dp

Dg2 = oo16pi2*( betag21 + oo16pi2 * betag22 )

Else
Dg2 = oo16pi2* betag21
End If

!--------------------
! g3
!--------------------

betag31  = -7._dp*(g3p3)

If (TwoLoopRGE) Then
betag32 = -(g3p3*(-11._dp*(g1p2) - 45._dp*(g2p2) + 260._dp*(g3p2) + 20._dp*(TrYdadjYd) +        &
&  20._dp*(TrYuadjYu)))/10._dp

betag33 = (65._dp*g3**7)/2._dp - 40*g3**5*Yt**2 + 15*g3**3*Yt**4

Dg3 = oo16pi2*( betag31 + oo16pi2 * betag32 + oo16pi2**2*betag33 )

Else
Dg3 = oo16pi2* betag31
End If

!--------------------
! Lam
!--------------------

betaLam1  = 27._dp*(g1p4)/100._dp + (9*g1p2*g2p2)/10._dp + 9._dp*(g2p4)               &
& /4._dp + 12._dp*(Lamp2) - 12._dp*(TrYdadjYdYdadjYd) - 4._dp*(TrYeadjYeYeadjYe)         &
&  - 12._dp*(TrYuadjYuYuadjYu) - (9*g1p2*Lam)/5._dp - 9*g2p2*Lam + 12*TrYdadjYd*Lam +    &
&  4*TrYeadjYe*Lam + 12*TrYuadjYu*Lam

If (TwoLoopRGE) Then
betaLam2 = -3411._dp*(g1p6)/1000._dp - (1677*g1p4*g2p2)/200._dp - (289*g1p2*g2p4)/40._dp +       &
&  305._dp*(g2p6)/8._dp + (54*g1p2*Lamp2)/5._dp + 54*g2p2*Lamp2 - 78._dp*(Lamp3) +       &
&  (8*g1p2*TrYdadjYdYdadjYd)/5._dp - 64*g3p2*TrYdadjYdYdadjYd + 60._dp*(TrYdadjYdYdadjYdYdadjYd) -&
&  24._dp*(TrYdadjYdYdadjYuYuadjYd) + 12._dp*(TrYdadjYuYuadjYdYdadjYd) - 12._dp*(TrYdadjYuYuadjYuYuadjYd) -&
&  (24*g1p2*TrYeadjYeYeadjYe)/5._dp + 20._dp*(TrYeadjYeYeadjYeYeadjYe) - (171*g1p4*TrYuadjYu)/50._dp +&
&  (63*g1p2*g2p2*TrYuadjYu)/5._dp - (9*g2p4*TrYuadjYu)/2._dp - 72*Lamp2*TrYuadjYu -      &
&  (16*g1p2*TrYuadjYuYuadjYu)/5._dp - 64*g3p2*TrYuadjYuYuadjYu + 60._dp*(TrYuadjYuYuadjYuYuadjYu) +&
&  (1887*g1p4*Lam)/200._dp + (117*g1p2*g2p2*Lam)/20._dp - (73*g2p4*Lam)/8._dp -          &
&  3*TrYdadjYdYdadjYd*Lam - 42*TrYdadjYuYuadjYd*Lam - TrYeadjYeYeadjYe*Lam +             &
&  (17*g1p2*TrYuadjYu*Lam)/2._dp + (45*g2p2*TrYuadjYu*Lam)/2._dp + 80*g3p2*TrYuadjYu*Lam -&
&  3*TrYuadjYuYuadjYu*Lam + (TrYdadjYd*(9._dp*(g1p4) - 45._dp*(g2p4) + 225*g2p2*Lam +    &
&  80*(10._dp*(g3p2) - 9._dp*(Lam))*Lam + g1p2*(54._dp*(g2p2) + 25._dp*(Lam))))/10._dp - &
&  (3*TrYeadjYe*(15._dp*(g1p4) - g1p2*(22._dp*(g2p2) + 25._dp*(Lam)) + 5*(g2p4 +         &
&  16._dp*(Lamp2) - 5*g2p2*Lam)))/10._dp

betaLam3 =   (5238*Lam**3*Yt**2 + Lam*Yt**2*(24*g3**2*Yt**2*(895 - 1296*Z3) + 27*Yt**4*(13 - 176*Z3) &
& + 32*g3**4*(311 - 36*Z3)) + 36*Lam**4*(299 + 168*Z3) + &
& 2*Yt**4*(16*g3**4*(-133 + 48*Z3) - 9*Yt**4*(533 + 96*Z3) + 48*g3**2*Yt**2*(-19 + 120*Z3)) + &
& 54*Lam**2*Yt**2*(16*g3**2*(-17 + 16*Z3) + Yt**2*(191 + 168*Z3)))/48._dp

DLam = oo16pi2*( betaLam1 + oo16pi2 * betaLam2 + oo16pi2**2*betaLam3 )

Else
DLam = oo16pi2* betaLam1
End If

Call Chop(DLam)

!--------------------
! Yu
!--------------------

betaYu1  = (-17._dp*(g1p2)/20._dp - 9._dp*(g2p2)/4._dp - 8._dp*(g3p2) +               &
&  3._dp*(TrYdadjYd) + TrYeadjYe + 3._dp*(TrYuadjYu))*Yu - (3*(YuadjYdYd -               &
&  YuadjYuYu))/2._dp

If (TwoLoopRGE) Then
betaYu2 = ((1187._dp*(g1p4) - 270*g1p2*g2p2 - 3450._dp*(g2p4) + 760*g1p2*g3p2 + 5400*g2p2*g3p2 -&
&  64800._dp*(g3p4) + 900._dp*(Lamp2) + 375*(g1p2 + 9._dp*(g2p2) + 32._dp*(g3p2))*TrYdadjYd -&
&  4050._dp*(TrYdadjYdYdadjYd) + 900._dp*(TrYdadjYuYuadjYd) + 1125*(g1p2 +               &
&  g2p2)*TrYeadjYe - 1350._dp*(TrYeadjYeYeadjYe) + 1275*g1p2*TrYuadjYu + 3375*g2p2*TrYuadjYu +&
&  12000*g3p2*TrYuadjYu - 4050._dp*(TrYuadjYuYuadjYu))*Yu)/600._dp + ((-43._dp*(g1p2) +  &
&  45._dp*(g2p2) - 1280._dp*(g3p2) + 300._dp*(TrYdadjYd) + 100._dp*(TrYeadjYe) +         &
&  300._dp*(TrYuadjYu))*YuadjYdYd + 20*(11._dp*(YuadjYdYdadjYdYd) - YuadjYdYdadjYuYu -   &
&  4._dp*(YuadjYuYuadjYdYd) + 6._dp*(YuadjYuYuadjYuYu)) + YuadjYuYu*(223._dp*(g1p2) +    &
&  675._dp*(g2p2) + 1280._dp*(g3p2) - 540._dp*(TrYdadjYd) - 180._dp*(TrYeadjYe) -        &
&  540._dp*(TrYuadjYu) - 480._dp*(Lam)))/80._dp

betaYu3 =  (Yt*(48*g3**2*(8*Lam*Yt**2 - 157*Yt**4) + 8*g3**4*Yt**2*(3827 - 1368*Z3) + &
& 32*g3**6*(-2083 + 960*Z3) + 9*(-24*Lam**3 + 5*Lam**2*Yt**2 + 528*Lam*Yt**4 + &
& 2*Yt**6*(113 + 36*Z3))))/48._dp

DYu = oo16pi2*( betaYu1 + oo16pi2 * betaYu2 )
DYu(3,3) = DYu(3,3) +  oo16pi2**3*betaYu3

Else
DYu = oo16pi2* betaYu1
End If

Call Chop(DYu)

!--------------------
! Yd
!--------------------

betaYd1  = (-((g1p2 + 9._dp*(g2p2) + 32._dp*(g3p2) - 12._dp*(TrYdadjYd)               &
&  - 4._dp*(TrYeadjYe) - 12._dp*(TrYuadjYu))*Yd) + 6*(YdadjYdYd - YdadjYuYu))/4._dp

If (TwoLoopRGE) Then
betaYd2 = (-127._dp*(g1p4)/600._dp - (27*g1p2*g2p2)/20._dp - 23._dp*(g2p4)/4._dp +              &
&  (31*g1p2*g3p2)/15._dp + 9*g2p2*g3p2 - 108._dp*(g3p4) + 3._dp*(Lamp2)/2._dp +          &
&  (5*(g1p2 + 9._dp*(g2p2) + 32._dp*(g3p2))*TrYdadjYd)/8._dp - 27._dp*(TrYdadjYdYdadjYd)/4._dp +&
&  3._dp*(TrYdadjYuYuadjYd)/2._dp + (15*(g1p2 + g2p2)*TrYeadjYe)/8._dp - 9._dp*(TrYeadjYeYeadjYe)/4._dp +&
&  (17*g1p2*TrYuadjYu)/8._dp + (45*g2p2*TrYuadjYu)/8._dp + 20*g3p2*TrYuadjYu -           &
&  27._dp*(TrYuadjYuYuadjYu)/4._dp)*Yd + ((-79._dp*(g1p2) + 45._dp*(g2p2) -              &
&  1280._dp*(g3p2) + 300._dp*(TrYdadjYd) + 100._dp*(TrYeadjYe) + 300._dp*(TrYuadjYu))*YdadjYuYu +&
&  20*(6._dp*(YdadjYdYdadjYdYd) - 4._dp*(YdadjYdYdadjYuYu) - YdadjYuYuadjYdYd +          &
&  11._dp*(YdadjYuYuadjYuYu)) + YdadjYdYd*(187._dp*(g1p2) + 675._dp*(g2p2) +             &
&  1280._dp*(g3p2) - 540._dp*(TrYdadjYd) - 180._dp*(TrYeadjYe) - 540._dp*(TrYuadjYu) -   &
&  480._dp*(Lam)))/80._dp

DYd = oo16pi2*( betaYd1 + oo16pi2 * betaYd2 )

Else
DYd = oo16pi2* betaYd1
End If

Call Chop(DYd)

!--------------------
! Ye
!--------------------

betaYe1  = (-9._dp*(g1p2)/4._dp - 9._dp*(g2p2)/4._dp + 3._dp*(TrYdadjYd)              &
&  + TrYeadjYe + 3._dp*(TrYuadjYu))*Ye + 3._dp*(YeadjYeYe)/2._dp

If (TwoLoopRGE) Then
betaYe2 = ((2742._dp*(g1p4) + 540*g1p2*g2p2 - 2300._dp*(g2p4) + 600._dp*(Lamp2) +               &
&  250*(g1p2 + 9._dp*(g2p2) + 32._dp*(g3p2))*TrYdadjYd - 2700._dp*(TrYdadjYdYdadjYd) +   &
&  600._dp*(TrYdadjYuYuadjYd) + 750*(g1p2 + g2p2)*TrYeadjYe - 900._dp*(TrYeadjYeYeadjYe) +&
&  850*g1p2*TrYuadjYu + 2250*g2p2*TrYuadjYu + 8000*g3p2*TrYuadjYu - 2700._dp*(TrYuadjYuYuadjYu))*Ye +&
&  15*(40._dp*(YeadjYeYeadjYeYe) + YeadjYeYe*(129._dp*(g1p2) + 225._dp*(g2p2) -          &
&  180._dp*(TrYdadjYd) - 60._dp*(TrYeadjYe) - 180._dp*(TrYuadjYu) - 160._dp*(Lam))))/400._dp

DYe = oo16pi2*( betaYe1 + oo16pi2 * betaYe2 )

Else
DYe = oo16pi2* betaYe1
End If

Call Chop(DYe)

!--------------------
! Mu
!--------------------

betaMu1  = (-9*g1p2*Mu)/10._dp - (9*g2p2*Mu)/2._dp + 6*Mu*TrYdadjYd + 2*Mu*TrYeadjYe +&
&  6*Mu*TrYuadjYu + 6*Mu*Lam

If (TwoLoopRGE) Then
betaMu2 = (Mu*(1671._dp*(g1p4) + 450*g1p2*g2p2 - 3625._dp*(g2p4) - 6000._dp*(Lamp2) -           &
&  5400._dp*(TrYdadjYdYdadjYd) - 8400._dp*(TrYdadjYuYuadjYd) - 1800._dp*(TrYeadjYeYeadjYe) +&
&  1700*g1p2*TrYuadjYu + 4500*g2p2*TrYuadjYu + 16000*g3p2*TrYuadjYu - 5400._dp*(TrYuadjYuYuadjYu) +&
&  100*TrYdadjYd*(5._dp*(g1p2) + 45._dp*(g2p2) + 160._dp*(g3p2) - 144._dp*(Lam)) +       &
&  300*TrYeadjYe*(5._dp*(g1p2) + 5._dp*(g2p2) - 16._dp*(Lam)) + 2880*g1p2*Lam +          &
&  14400*g2p2*Lam - 14400*TrYuadjYu*Lam))/400._dp

DMu = oo16pi2*( betaMu1 + oo16pi2 * betaMu2 )

Else
DMu = oo16pi2* betaMu1
End If

Call Chop(DMu)

!--------------------
! v
!--------------------

betav1  = (v*(9._dp*(g1p2) + 45._dp*(g2p2) - 60._dp*(TrYdadjYd) - 20._dp*(TrYeadjYe)  &
&  - 60._dp*(TrYuadjYu) + 3*g1p2*Xi + 15*g2p2*Xi))/20._dp

If (TwoLoopRGE) Then
betav2 = (v*(-1293._dp*(g1p4) - 270*g1p2*g2p2 + 6775._dp*(g2p4) - 1200._dp*(Lamp2) +           &
&  5400._dp*(TrYdadjYdYdadjYd) - 1200._dp*(TrYdadjYuYuadjYd) + 1800._dp*(TrYeadjYeYeadjYe) -&
&  1700*g1p2*TrYuadjYu - 4500*g2p2*TrYuadjYu - 16000*g3p2*TrYuadjYu + 5400._dp*(TrYuadjYuYuadjYu) +&
&  18*g1p4*Xi + 180*g1p2*g2p2*Xi + 2250*g2p4*Xi - 360*g1p2*TrYuadjYu*Xi - 1800*g2p2*TrYuadjYu*Xi -&
&  60*TrYeadjYe*(5*g2p2*(5 + 2._dp*(Xi)) + g1p2*(25 + 2._dp*(Xi))) - 20*TrYdadjYd*(800._dp*(g3p2) +&
&  45*g2p2*(5 + 2._dp*(Xi)) + g1p2*(25 + 18._dp*(Xi))) + 18*g1p4*Xip2 + 180*g1p2*g2p2*Xip2 -&
&  450*g2p4*Xip2))/800._dp

Dv = oo16pi2*( betav1 + oo16pi2 * betav2 )

Else
Dv = oo16pi2* betav1
End If

Call ParametersToG62_SM(Dg1,Dg2,Dg3,DLam,DYu,DYd,DYe,DMu,Dv,f)

Iname = Iname - 1

End Subroutine rge62_SM

End Module RGEs_SM_HC


text:
subroutine other(res, var, nz, ny, nx)

    use kind

    implicit none

    integer(kind=ni) :: nx, ny, nz

    real   (kind=nr), intent(in)  :: var (0:nx, 0:ny, 0:nz)
    real   (kind=nr), intent(out) :: res (0:nx, 0:ny, 0:nz)

    integer(kind=ni) :: i, j, k

    do k = 0, nz
      do j = 0, ny
        do i = 0, nx
          res(i, j, k) = var(i, j, k) * var(i, j, k)
        enddo
      enddo
    enddo
end subroutine other

max_stars_repo_name: dennissergeev/dummyforpy

max_issues_repo_name: dennissergeev/dummyforpy

max_forks_repo_name: dennissergeev/dummyforpy


text:
      subroutine triangle7new(k1,k2,k3,k4,k5,k6,za,zb,app,apm)
      implicit none
      include 'constants.f'
      include 'masses.f'
      include 'sprods_com.f'
      include 'zprods_decl.f'
      include 'Higgsint.f'
      double precision mtsq
      integer k1,k2,k3,k4,k5,k6,i1,i2
      double complex app,apm,triamp,za34b,iza34b,iza,izb,
     & t2sum,t2prod,t1234sum
c--- statement functions
      iza(i1,i2)=cone/za(i1,i2)
      izb(i1,i2)=cone/zb(i1,i2)
      za34b(i1,i2)=za(i1,k3)*zb(k3,i2)+za(i1,k4)*zb(k4,i2)
      iza34b(i1,i2)=cone/za34b(i1,i2)
c--- end statement functions

      mtsq=mt**2
      apm=czip

c--- Triangle 7
      t2sum=-(-za(k1,k2)*zb(k2,k1)/za34b(k1,k1)*za(k3,k4)*zb(k4,k3)
     &        -za(k1,k2)*zb(k2,k1)/za34b(k1,k1)*mtsq
     &        +za34b(k2,k2))
      t2prod=-za(k1,k2)*zb(k2,k1)/za34b(k1,k1)*mtsq
     & *(za34b(k2,k2)
     &  -za(k3,k4)*zb(k4,k3)/za34b(k1,k1)*za(k1,k2)*zb(k2,k1))
      t1234sum=-(za(k1,k2)*zb(k2,k1)/za34b(k1,k1)*za(k3,k4)*zb(k4,k3)
     &          +za(k1,k2)*zb(k2,k1)/za34b(k1,k1)*mtsq
     &          +za(k1,k2)*zb(k2,k1))
c      t1234prod=-za(k1,k2)*zb(k2,k1)/za34b(k1,k1)*mtsq
c     & *(za34b(k2,k2)
c     &  -za(k3,k4)*zb(k4,k3)/za34b(k1,k1)*za(k1,k2)*zb(k2,k1))
      triamp =  + mtsq * (  - za(k1,k3)**2*za(k1,k5)*zb(k1,k2)*zb(k2,k6
     &    )*zb(k3,k4)*iza(k1,k2)**2*za34b(k2,k1)*iza34b(k1,k1)*iza34b(
     &    k1,k2) - za(k1,k3)**2*za(k1,k5)*zb(k1,k2)*zb(k2,k6)*zb(k3,k4)
     &    *iza(k1,k2)**2*iza34b(k1,k2)**2*t2sum - za(k1,k3)**2*za(k1,k5
     &    )*zb(k1,k2)*zb(k3,k4)*iza(k1,k2)**3*za34b(k1,k6)*iza34b(k1,k2
     &    )**2*t1234sum - za(k1,k3)**2*za(k1,k5)*zb(k1,k2)*zb(k3,k4)*
     &    iza(k1,k2)**3*za34b(k1,k6)*iza34b(k1,k2)**2*t2sum - za(k1,k3)
     &    **2*za(k1,k5)*zb(k1,k6)*zb(k3,k4)*iza(k1,k2)**3*za34b(k2,k1)*
     &    iza34b(k1,k1) - za(k1,k3)**2*za(k1,k5)*zb(k1,k6)*zb(k3,k4)*
     &    iza(k1,k2)**3*iza34b(k1,k2)*t2sum + za(k1,k3)*za(k1,k5)*za(k2
     &    ,k3)*zb(k1,k2)*zb(k2,k6)*zb(k3,k4)*iza(k1,k2)**2*iza34b(k1,k2
     &    ) + za(k1,k3)*za(k1,k5)*za(k2,k3)*zb(k1,k6)*zb(k3,k4)*iza(k1,
     &    k2)**3 + za(k1,k3)*za(k1,k5)*za(k3,k4)*zb(k1,k2)*zb(k1,k4)*
     &    zb(k2,k6)*zb(k3,k4)*iza(k1,k2)*iza34b(k1,k1)*iza34b(k1,k2) +
     &    za(k1,k3)*za(k1,k5)*za(k3,k4)*zb(k1,k4)*zb(k1,k6)*zb(k3,k4)*
     &    iza(k1,k2)**2*iza34b(k1,k1) )
      triamp = triamp + mtsq * (  - za(k1,k3)*za(k1,k5)*zb(k1,k4)*zb(k2
     &    ,k6)*iza(k1,k2)**2*za34b(k2,k1)*iza34b(k1,k1) - za(k1,k3)*za(
     &    k1,k5)*zb(k1,k4)*zb(k2,k6)*iza(k1,k2)**2*iza34b(k1,k2)*t2sum
     &     - 2.D0*za(k1,k3)*za(k1,k5)*zb(k1,k4)*iza(k1,k2)**3*za34b(k1,
     &    k6)*za34b(k2,k1)*iza34b(k1,k1) + za(k1,k3)*za(k1,k5)*zb(k1,k4
     &    )*iza(k1,k2)**3*za34b(k1,k6)*za34b(k2,k2)*iza34b(k1,k2) - za(
     &    k1,k3)*za(k1,k5)*zb(k1,k4)*iza(k1,k2)**3*za34b(k1,k6)*iza34b(
     &    k1,k2)*t2sum + za(k1,k3)*za(k2,k5)*zb(k1,k4)*zb(k2,k6)*iza(k1
     &    ,k2)**2 + za(k1,k3)*za(k2,k5)*zb(k1,k4)*iza(k1,k2)**3*za34b(
     &    k1,k6) + za(k1,k3)*zb(k1,k4)*iza(k1,k2)**2*za34b(k1,k6)*
     &    za34b(k5,k1)*iza34b(k1,k1) - za(k1,k3)*zb(k1,k4)*iza(k1,k2)**
     &    2*za34b(k1,k6)*za34b(k5,k2)*iza34b(k1,k2) + 2.D0*za(k1,k5)*
     &    za(k3,k4)*zb(k1,k4)**2*zb(k1,k6)*iza(k1,k2)**2*izb(k1,k2)*
     &    za34b(k1,k2)*iza34b(k1,k1) - za(k1,k5)*za(k3,k4)*zb(k1,k4)**2
     &    *zb(k1,k6)*iza(k1,k2)**2*izb(k1,k2)*za34b(k2,k2)*iza34b(k2,k1
     &    ) )
      triamp = triamp + mtsq * ( za(k1,k5)*za(k3,k4)*zb(k1,k4)**2*zb(k1
     &    ,k6)*iza(k1,k2)**2*izb(k1,k2)*iza34b(k2,k1)*t2sum + 2.D0*za(
     &    k1,k5)*za(k3,k4)*zb(k1,k4)**2*zb(k2,k6)*iza(k1,k2)*iza34b(k1,
     &    k1) + 2.D0*za(k1,k5)*za(k3,k4)*zb(k1,k4)**2*iza(k1,k2)**2*
     &    za34b(k1,k6)*iza34b(k1,k1) - za(k2,k5)*za(k3,k4)*zb(k1,k4)**2
     &    *zb(k1,k6)*iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k1)*iza34b(k2,k1
     &    )**2*t2sum - za(k2,k5)*za(k3,k4)*zb(k1,k4)**2*zb(k1,k6)*iza(
     &    k1,k2)**2*izb(k1,k2)*za34b(k1,k2)*iza34b(k2,k1) - za(k2,k5)*
     &    za(k3,k4)*zb(k1,k4)**2*zb(k2,k6)*iza(k1,k2)*iza34b(k2,k1) -
     &    za(k2,k5)*za(k3,k4)*zb(k1,k4)**2*iza(k1,k2)**2*za34b(k1,k6)*
     &    iza34b(k2,k1) - za(k3,k4)*zb(k1,k4)**2*zb(k1,k6)*iza(k1,k2)*
     &    izb(k1,k2)*za34b(k1,k2)*za34b(k5,k1)*iza34b(k1,k1)*iza34b(k2,
     &    k1) + za(k3,k4)*zb(k1,k4)**2*zb(k1,k6)*iza(k1,k2)*izb(k1,k2)*
     &    za34b(k5,k1)*iza34b(k2,k1)**2*t1234sum + za(k3,k4)*zb(k1,k4)
     &    **2*zb(k1,k6)*iza(k1,k2)*izb(k1,k2)*za34b(k5,k2)*iza34b(k2,k1
     &    ) )
      triamp = triamp + mtsq * (  - za(k3,k4)*zb(k1,k4)**2*iza(k1,k2)*
     &    za34b(k1,k6)*za34b(k5,k1)*iza34b(k1,k1)*iza34b(k2,k1) )
      triamp = triamp + mtsq**2 * (  - za(k1,k3)*za(k1,k5)*zb(k1,k2)*
     &    zb(k1,k4)*zb(k2,k6)*iza(k1,k2)*iza34b(k1,k1)*iza34b(k1,k2) -
     &    2.D0*za(k1,k3)*za(k1,k5)*zb(k1,k4)*zb(k1,k6)*iza(k1,k2)**2*
     &    iza34b(k1,k1) + za(k1,k3)*za(k2,k5)*zb(k1,k4)*zb(k1,k6)*iza(
     &    k1,k2)**2*iza34b(k2,k1) + za(k1,k3)*zb(k1,k4)*zb(k1,k6)*iza(
     &    k1,k2)*za34b(k5,k1)*iza34b(k1,k1)*iza34b(k2,k1) )
      triamp = triamp + za(k1,k3)**2*za(k1,k5)*zb(k2,k6)*zb(k3,k4)*iza(
     & k1,k2)**3*za34b(k1,k1)*iza34b(k1,k2)**2*t2prod - za(k1,k3)**2*
     &    za(k1,k5)*zb(k2,k6)*zb(k3,k4)*iza(k1,k2)**3*za34b(k1,k1)*
     &    iza34b(k1,k2)**2*t2sum**2 - za(k1,k3)**2*za(k1,k5)*zb(k2,k6)*
     &    zb(k3,k4)*iza(k1,k2)**3*za34b(k2,k1)*iza34b(k1,k2)*t2sum +
     &    za(k1,k3)**2*za(k1,k5)*zb(k3,k4)*iza(k1,k2)**4*za34b(k1,k1)*
     &    za34b(k1,k6)*iza34b(k1,k2)**2*t2prod - za(k1,k3)**2*za(k1,k5)
     &    *zb(k3,k4)*iza(k1,k2)**4*za34b(k1,k1)*za34b(k1,k6)*iza34b(k1,
     &    k2)**2*t2sum**2 - za(k1,k3)**2*za(k1,k5)*zb(k3,k4)*iza(k1,k2)
     &    **4*za34b(k1,k6)*za34b(k2,k1)*iza34b(k1,k2)*t2sum + za(k1,k3)
     &    *za(k1,k5)*za(k2,k3)*zb(k2,k6)*zb(k3,k4)*iza(k1,k2)**3*za34b(
     &    k1,k1)*iza34b(k1,k2)*t2sum + za(k1,k3)*za(k1,k5)*za(k2,k3)*
     &    zb(k3,k4)*iza(k1,k2)**4*za34b(k1,k1)*za34b(k1,k6)*iza34b(k1,
     &    k2)*t2sum + za(k1,k3)*za(k1,k5)*za(k3,k4)*zb(k1,k4)*zb(k2,k6)
     &    *zb(k3,k4)*iza(k1,k2)**2*za34b(k2,k1)*iza34b(k1,k1) + za(k1,
     &    k3)*za(k1,k5)*za(k3,k4)*zb(k1,k4)*zb(k2,k6)*zb(k3,k4)*iza(k1,
     &    k2)**2*iza34b(k1,k2)*t2sum
      triamp = triamp + za(k1,k3)*za(k1,k5)*za(k3,k4)*zb(k1,k4)*zb(k3,
     & k4)*iza(k1,k2)**3*za34b(k1,k6)*za34b(k2,k1)*iza34b(k1,k1) + za(
     &    k1,k3)*za(k1,k5)*za(k3,k4)*zb(k1,k4)*zb(k3,k4)*iza(k1,k2)**3*
     &    za34b(k1,k6)*iza34b(k1,k2)*t2sum - za(k1,k5)*za(k2,k3)*za(k3,
     &    k4)*zb(k1,k4)*zb(k2,k6)*zb(k3,k4)*iza(k1,k2)**2 - za(k1,k5)*
     &    za(k2,k3)*za(k3,k4)*zb(k1,k4)*zb(k3,k4)*iza(k1,k2)**3*za34b(
     &    k1,k6)

      app=triamp*(-half*im)/(s(k3,k4)*s(k5,k6))
      if (Higgsint) return
      triamp =  + mtsq * (  - za(k1,k3)**2*za(k1,k5)*zb(k1,k6)*zb(k3,k4
     &    )*iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k1)*iza34b(k1,k2)**3*
     &    t1234sum**2 + za(k1,k3)**2*za(k1,k5)*zb(k1,k6)*zb(k3,k4)*iza(
     &    k1,k2)**2*izb(k1,k2)*za34b(k1,k1)*iza34b(k1,k2)**3*t2sum**2
     &     + za(k1,k3)**2*za(k1,k5)*zb(k1,k6)*zb(k3,k4)*iza(k1,k2)**2*
     &    izb(k1,k2)*za34b(k2,k1)**2*iza34b(k1,k1)*iza34b(k1,k2) + za(
     &    k1,k3)**2*za(k1,k5)*zb(k1,k6)*zb(k3,k4)*iza(k1,k2)**2*izb(k1,
     &    k2)*za34b(k2,k1)*iza34b(k1,k2)**2*t1234sum + 2.D0*za(k1,k3)**
     &    2*za(k1,k5)*zb(k1,k6)*zb(k3,k4)*iza(k1,k2)**2*izb(k1,k2)*
     &    za34b(k2,k1)*iza34b(k1,k2)**2*t2sum + za(k1,k3)**2*za(k2,k5)*
     &    zb(k1,k6)*zb(k3,k4)*iza(k1,k2)*za34b(k2,k1)*iza34b(k1,k1)*
     &    iza34b(k1,k2) - za(k1,k3)**2*za(k2,k5)*zb(k1,k6)*zb(k3,k4)*
     &    iza(k1,k2)*iza34b(k1,k2)**2*t1234sum - za(k1,k3)*za(k1,k5)*
     &    za(k2,k3)*zb(k1,k6)*zb(k3,k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(
     &    k1,k1)*iza34b(k1,k2)**2*t1234sum - za(k1,k3)*za(k1,k5)*za(k2,
     &    k3)*zb(k1,k6)*zb(k3,k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k1)
     &    *iza34b(k1,k2)**2*t2sum )
      triamp = triamp + mtsq * (  - za(k1,k3)*za(k1,k5)*za(k2,k3)*zb(k1
     &    ,k6)*zb(k3,k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(k2,k1)*iza34b(
     &    k1,k2) - za(k1,k3)*za(k1,k5)*za(k3,k4)*zb(k1,k4)*zb(k1,k6)*
     &    zb(k3,k4)*iza(k1,k2)*izb(k1,k2)*za34b(k2,k1)*iza34b(k1,k1)*
     &    iza34b(k1,k2) - za(k1,k3)*za(k1,k5)*za(k3,k4)*zb(k1,k4)*zb(k1
     &    ,k6)*zb(k3,k4)*iza(k1,k2)*izb(k1,k2)*iza34b(k1,k2)**2*
     &    t1234sum - za(k1,k3)*za(k1,k5)*za(k3,k4)*zb(k1,k4)*zb(k1,k6)*
     &    zb(k3,k4)*iza(k1,k2)*izb(k1,k2)*iza34b(k1,k2)**2*t2sum - za(
     &    k1,k3)*za(k1,k5)*zb(k1,k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(k1,
     &    k1)*za34b(k1,k6)*iza34b(k1,k2)**3*t1234sum**2 + za(k1,k3)*za(
     &    k1,k5)*zb(k1,k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k1)*za34b(
     &    k1,k6)*iza34b(k1,k2)**3*t2sum**2 + za(k1,k3)*za(k1,k5)*zb(k1,
     &    k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k6)*za34b(k2,k1)**2*
     &    iza34b(k1,k1)*iza34b(k1,k2) + za(k1,k3)*za(k1,k5)*zb(k1,k4)*
     &    iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k6)*za34b(k2,k1)*iza34b(k1,
     &    k2)**2*t1234sum )
      triamp = triamp + mtsq * ( 2.D0*za(k1,k3)*za(k1,k5)*zb(k1,k4)*
     &    iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k6)*za34b(k2,k1)*iza34b(k1,
     &    k2)**2*t2sum - za(k1,k3)*za(k2,k3)*za(k2,k5)*zb(k1,k6)*zb(k3,
     &    k4)*iza(k1,k2)*iza34b(k1,k2) - za(k1,k3)*za(k2,k5)*za(k3,k4)*
     &    zb(k1,k4)*zb(k1,k6)*zb(k3,k4)*iza34b(k1,k1)*iza34b(k1,k2) -
     &    za(k1,k3)*za(k2,k5)*zb(k1,k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(
     &    k1,k1)*za34b(k1,k6)*iza34b(k1,k2)**2*t1234sum - za(k1,k3)*za(
     &    k2,k5)*zb(k1,k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k1)*za34b(
     &    k1,k6)*iza34b(k1,k2)**2*t2sum - za(k1,k3)*za(k2,k5)*zb(k1,k4)
     &    *iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k6)*za34b(k2,k1)*iza34b(k1
     &    ,k2) + za(k1,k3)*za(k2,k5)*zb(k1,k4)*iza(k1,k2)*za34b(k1,k6)*
     &    za34b(k2,k1)*iza34b(k1,k1)*iza34b(k1,k2) - za(k1,k3)*za(k2,k5
     &    )*zb(k1,k4)*iza(k1,k2)*za34b(k1,k6)*iza34b(k1,k2)**2*t1234sum
     &     - za(k1,k5)*za(k3,k4)*zb(k1,k4)**2*zb(k1,k6)*iza(k1,k2)*izb(
     &    k1,k2)**2*za34b(k2,k1)*iza34b(k1,k1) - za(k1,k5)*za(k3,k4)*
     &    zb(k1,k4)**2*iza(k1,k2)*izb(k1,k2)*za34b(k1,k6)*za34b(k2,k1)*
     &    iza34b(k1,k1)*iza34b(k1,k2) )
      triamp = triamp + mtsq * (  - za(k1,k5)*za(k3,k4)*zb(k1,k4)**2*
     &    iza(k1,k2)*izb(k1,k2)*za34b(k1,k6)*iza34b(k1,k2)**2*t1234sum
     &     - za(k1,k5)*za(k3,k4)*zb(k1,k4)**2*iza(k1,k2)*izb(k1,k2)*
     &    za34b(k1,k6)*iza34b(k1,k2)**2*t2sum + za(k2,k5)*za(k3,k4)*zb(
     &    k1,k4)**2*zb(k1,k6)*iza(k1,k2)*izb(k1,k2)**2 - za(k2,k5)*za(
     &    k3,k4)*zb(k1,k4)**2*zb(k1,k6)*izb(k1,k2)*iza34b(k1,k1) - za(
     &    k2,k5)*za(k3,k4)*zb(k1,k4)**2*za34b(k1,k6)*iza34b(k1,k1)*
     &    iza34b(k1,k2) )
      triamp = triamp + mtsq**2 * ( za(k1,k3)*za(k1,k5)*zb(k1,k4)*zb(k1
     &    ,k6)*iza(k1,k2)*izb(k1,k2)*za34b(k2,k1)*iza34b(k1,k1)*iza34b(
     &    k1,k2) + za(k1,k3)*za(k1,k5)*zb(k1,k4)*zb(k1,k6)*iza(k1,k2)*
     &    izb(k1,k2)*iza34b(k1,k2)**2*t1234sum + za(k1,k3)*za(k1,k5)*
     &    zb(k1,k4)*zb(k1,k6)*iza(k1,k2)*izb(k1,k2)*iza34b(k1,k2)**2*
     &    t2sum + za(k1,k3)*za(k2,k5)*zb(k1,k4)*zb(k1,k6)*iza34b(k1,k1)
     &    *iza34b(k1,k2) )
      triamp = triamp - 2.D0*za(k1,k3)**2*za(k1,k5)*zb(k3,k4)*iza(k1,k2
     & )**3*izb(k1,k2)*za34b(k1,k1)**2*za34b(k1,k6)*iza34b(k1,k2)**4*
     & t1234sum*t2prod + za(k1,k3)**2*za(k1,k5)*zb(k3,k4)*iza(k1,k2)**3
     &    *izb(k1,k2)*za34b(k1,k1)**2*za34b(k1,k6)*iza34b(k1,k2)**4*
     &    t1234sum**3 - 2.D0*za(k1,k3)**2*za(k1,k5)*zb(k3,k4)*iza(k1,k2
     &    )**3*izb(k1,k2)*za34b(k1,k1)**2*za34b(k1,k6)*iza34b(k1,k2)**4
     &    *t2sum*t2prod + za(k1,k3)**2*za(k1,k5)*zb(k3,k4)*iza(k1,k2)**
     &    3*izb(k1,k2)*za34b(k1,k1)**2*za34b(k1,k6)*iza34b(k1,k2)**4*
     &    t2sum**3 - za(k1,k3)**2*za(k1,k5)*zb(k3,k4)*iza(k1,k2)**3*
     &    izb(k1,k2)*za34b(k1,k1)*za34b(k1,k6)*za34b(k2,k1)*iza34b(k1,
     &    k2)**3*t2prod - za(k1,k3)**2*za(k1,k5)*zb(k3,k4)*iza(k1,k2)**
     &    3*izb(k1,k2)*za34b(k1,k1)*za34b(k1,k6)*za34b(k2,k1)*iza34b(k1
     &    ,k2)**3*t1234sum**2 + 2.D0*za(k1,k3)**2*za(k1,k5)*zb(k3,k4)*
     &    iza(k1,k2)**3*izb(k1,k2)*za34b(k1,k1)*za34b(k1,k6)*za34b(k2,
     &    k1)*iza34b(k1,k2)**3*t2sum**2 + za(k1,k3)**2*za(k1,k5)*zb(k3,
     &    k4)*iza(k1,k2)**3*izb(k1,k2)*za34b(k1,k6)*za34b(k2,k1)**2*
     &    iza34b(k1,k2)**2*t2sum
      triamp = triamp - za(k1,k3)**2*za(k2,k5)*zb(k3,k4)*iza(k1,k2)**2*
     & za34b(k1,k1)*za34b(k1,k6)*iza34b(k1,k2)**3*t2prod + za(k1,k3)**2
     &    *za(k2,k5)*zb(k3,k4)*iza(k1,k2)**2*za34b(k1,k1)*za34b(k1,k6)*
     &    iza34b(k1,k2)**3*t1234sum**2 - za(k1,k3)**2*za(k2,k5)*zb(k3,
     &    k4)*iza(k1,k2)**2*za34b(k1,k6)*za34b(k2,k1)*iza34b(k1,k2)**2*
     &    t1234sum + za(k1,k3)*za(k1,k5)*za(k2,k3)*zb(k3,k4)*iza(k1,k2)
     &    **3*izb(k1,k2)*za34b(k1,k1)**2*za34b(k1,k6)*iza34b(k1,k2)**3*
     &    t1234sum**2 - za(k1,k3)*za(k1,k5)*za(k2,k3)*zb(k3,k4)*iza(k1,
     &    k2)**3*izb(k1,k2)*za34b(k1,k1)**2*za34b(k1,k6)*iza34b(k1,k2)
     &    **3*t2sum**2 - za(k1,k3)*za(k1,k5)*za(k2,k3)*zb(k3,k4)*iza(k1
     &    ,k2)**3*izb(k1,k2)*za34b(k1,k1)*za34b(k1,k6)*za34b(k2,k1)*
     &    iza34b(k1,k2)**2*t2sum + za(k1,k3)*za(k1,k5)*za(k3,k4)*zb(k1,
     &    k4)*zb(k3,k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k1)*za34b(k1,
     &    k6)*iza34b(k1,k2)**3*t1234sum**2 - za(k1,k3)*za(k1,k5)*za(k3,
     &    k4)*zb(k1,k4)*zb(k3,k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k1)
     &    *za34b(k1,k6)*iza34b(k1,k2)**3*t2sum**2
      triamp = triamp - za(k1,k3)*za(k1,k5)*za(k3,k4)*zb(k1,k4)*zb(k3,
     & k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k6)*za34b(k2,k1)**2*
     & iza34b(k1,k1)*iza34b(k1,k2) - za(k1,k3)*za(k1,k5)*za(k3,k4)*zb(
     &    k1,k4)*zb(k3,k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k6)*za34b(
     &    k2,k1)*iza34b(k1,k2)**2*t1234sum - 2.D0*za(k1,k3)*za(k1,k5)*
     &    za(k3,k4)*zb(k1,k4)*zb(k3,k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(
     &    k1,k6)*za34b(k2,k1)*iza34b(k1,k2)**2*t2sum + za(k1,k3)*za(k2,
     &    k3)*za(k2,k5)*zb(k3,k4)*iza(k1,k2)**2*za34b(k1,k1)*za34b(k1,
     &    k6)*iza34b(k1,k2)**2*t1234sum - za(k1,k3)*za(k2,k5)*za(k3,k4)
     &    *zb(k1,k4)*zb(k3,k4)*iza(k1,k2)*za34b(k1,k6)*za34b(k2,k1)*
     &    iza34b(k1,k1)*iza34b(k1,k2) + za(k1,k3)*za(k2,k5)*za(k3,k4)*
     &    zb(k1,k4)*zb(k3,k4)*iza(k1,k2)*za34b(k1,k6)*iza34b(k1,k2)**2*
     &    t1234sum + za(k1,k5)*za(k2,k3)*za(k3,k4)*zb(k1,k4)*zb(k3,k4)*
     &    iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k1)*za34b(k1,k6)*iza34b(k1,
     &    k2)**2*t1234sum + za(k1,k5)*za(k2,k3)*za(k3,k4)*zb(k1,k4)*zb(
     &    k3,k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k1)*za34b(k1,k6)*
     &    iza34b(k1,k2)**2*t2sum
      triamp = triamp + za(k1,k5)*za(k2,k3)*za(k3,k4)*zb(k1,k4)*zb(k3,
     & k4)*iza(k1,k2)**2*izb(k1,k2)*za34b(k1,k6)*za34b(k2,k1)*iza34b(k1
     & ,k2) + za(k2,k3)*za(k2,k5)*za(k3,k4)*zb(k1,k4)*zb(k3,k4)*iza(k1,
     &    k2)*za34b(k1,k6)*iza34b(k1,k2)

      apm=triamp*(-half*im)/(s(k3,k4)*s(k5,k6))
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
!
!*******************************************************************************
!
subroutine Legendre_expand(jx_max,xval,Ang_Dist,L_max,Ang_L)
!
!*******************************************************************************
!
!  Discussion:
!
!  FORTRAN 90 subroutine to find the coefficients for an expaansion into
!  Legendre polynomials
!
!  Ang_Dist(x) = Sum Ang_L(L)*poly(L,1,alpha, beta,x)
!
!  Ang_L(L) = int(-1,1) Ang_Dist(x)*poly(L,1,alpha,beta,x)*dx*(2L+1)/2.
!
!  Note, input array is interpolated to be put onto the grid for the
!  Gauss-Legendre integration points
!
!   Dependencies:
!
!     Modules:
!
!        None
!
!     Subroutines:
!
!        None
!
!     External functions:
!
!        real(kind=8) :: Poly
!        real(kind=8) :: interp
!
!     MPI routines:
!
!        None
!
!  Licensing:
!
!    SPDX-License-Identifier: MIT
!
!  Date:
!
!    11 May 2021
!
!  Author:
!
!      Erich Ormand, LLNL
!
!*******************************************************************************
!
   implicit none
   integer(kind=4), intent(in) :: jx_max
   real(kind=8), intent(in) :: xval(jx_max),Ang_Dist(jx_max)
   integer(kind=4), intent(in) :: L_max
   real(kind=8), intent(out) :: Ang_L(0:L_max)
!------------------------------------------------------
   integer(kind=4) :: L, LL_max
   real(kind=8) :: sum, value
   real(kind=8) :: alpha, beta
   real(kind=8) :: xnorm
!   integer(kind=4) :: jx
   integer(kind=4) :: n
   integer(kind=4) :: n_gleg
   real(kind=8), allocatable :: x_gleg(:), w_gleg(:)
!-------------   External function calls
   real(kind=8) :: Poly
   real(kind=8) :: interp
!------------------------------------------------------
!
!----   Set up data for Legendre integration
!----   Can't have more polynomials than number of data points in Ang_Dist
!----   The array Ang_Dist is assumed to be on a fixed grid
!----   with values xval(jx) and fixed spacing dx with dx = xval(2) - xval(1)
!----   First number of computed coefficients is LL_max
!----   set LL_max = L_max, but if L_max > jx_max, adjust
!----   Then the same for number of Gauss-Legendre points. Must be <= jx_max
!----   In general, preferably less
!
   LL_max = L_max
   if(L_max > jx_max)LL_max = jx_max - 1
   n_gleg = jx_max

   Ang_L(0) = 0.5d0
   if(LL_max == 0)return    !   All done, just one value and must be 0.5

   if(.not.allocated(x_gleg))allocate(x_gleg(n_gleg))
   if(.not.allocated(w_gleg))allocate(w_gleg(n_gleg))

   alpha = 0.0d0
   beta = 0.0d0
   call gauss_quad(n_gleg, 1, alpha, beta, x_gleg, w_gleg)

   xnorm = 0.0d0
   do n = 1, n_gleg
      value = interp(x_gleg(n),jx_max,xval,Ang_Dist)
      xnorm = xnorm + value*w_gleg(n)
   end do

   Ang_L(0:LL_max) = 0.0d0

   do L = 0, L_max
      sum = 0.0d0
      do n = 1, n_gleg
         value = interp(x_gleg(n),jx_max,xval,Ang_Dist)
         sum = sum + value*w_gleg(n)*poly(L,1,alpha,beta,x_gleg(n))
      end do
      Ang_L(L) = sum*0.5d0*(2.0d0*real(L,kind=8)+1.0d0)/xnorm
   end do
   if(allocated(x_gleg))deallocate(x_gleg)
   if(allocated(w_gleg))deallocate(w_gleg)
   return
end subroutine Legendre_expand

max_stars_repo_name: LLNL/Yet-Another-Hauser-Feshbach-Code

max_issues_repo_name: LLNL/Yet-Another-Hauser-Feshbach-Code

max_forks_repo_name: LLNL/Yet-Another-Hauser-Feshbach-Code