text:
!----------------------------------------------------------------
! (c) Copyright, 2016 by the Regents of the University of California.
! Quadrupoleclass: Quadrupole beam line element class
!             in Lattice module of APPLICATION layer.
!
! MODULE  : ... Quadrupoleclass
! VERSION : ... 1.0
!> @author
!> Ji Qiang, LBNL
!> @date
!> 2016
!
! DESCRIPTION:
!> This class defines the linear transfer map and field
!> for the quadrupole beam line elment.
! Comments:
!----------------------------------------------------------------
      module Quadrupoleclass
        use PhysConstclass
        use Dataclass
        integer, private, parameter :: Nparam = 9
        type Quadrupole
          !Itype = 1
          integer :: Nseg,Mapstp,Itype !< Itype = 1
          double precision :: Length
          double precision, dimension(Nparam) :: Param
          ! Param(1) : zedge
          !      (2) : quad gradient
          !      (3) : file ID
          !      (4) : radius
          !      (5) : x misalignment error
          !      (6) : y misalignment error
          !      (7) : rotation error x
          !      (8) : rotation error y
          !      (9) : rotation error z
        end type Quadrupole
        interface getparam_Quadrupole
          module procedure getparam1_Quadrupole,  &
                          getparam2_Quadrupole,   &
                          getparam3_Quadrupole
        end interface
        interface setparam_Quadrupole
          module procedure setparam1_Quadrupole,  &
                          setparam2_Quadrupole, setparam3_Quadrupole
        end interface
      contains
        subroutine construct_Quadrupole(this,numseg,nmpstp,type,blength)
        implicit none
        type (Quadrupole), intent(out) :: this
        integer, intent(in) :: numseg,nmpstp,type
        double precision, intent(in) :: blength

        this%Nseg = numseg
        this%Mapstp = nmpstp
        this%Itype = type
        this%Length = blength

        this%Param = 0.0

        end subroutine construct_Quadrupole

        subroutine setparam1_Quadrupole(this,i,value)
        implicit none
        type (Quadrupole), intent(out) :: this
        integer, intent(in) :: i
        double precision, intent(in) :: value

        this%Param(i) = value

        end subroutine setparam1_Quadrupole

        subroutine setparam2_Quadrupole(this,values)
        implicit none
        type (Quadrupole), intent(out) :: this
        double precision, dimension(:), intent(in) :: values

        this%Param = values

        end subroutine setparam2_Quadrupole

        subroutine setparam3_Quadrupole(this,numseg,nmpstp,type,blength)
        implicit none
        type (Quadrupole), intent(out) :: this
        integer, intent(in) :: numseg,nmpstp,type
        double precision, intent(in) :: blength

        this%Nseg = numseg
        this%Mapstp = nmpstp
        this%Itype = type
        this%Length = blength

        end subroutine setparam3_Quadrupole

        subroutine getparam1_Quadrupole(this,i,blparam)
        implicit none
        type (Quadrupole), intent(in) :: this
        integer, intent(in) :: i
        double precision, intent(out) :: blparam

        blparam = this%Param(i)

        end subroutine getparam1_Quadrupole

        subroutine getparam2_Quadrupole(this,blparams)
        implicit none
        type (Quadrupole), intent(in) :: this
        double precision, dimension(:), intent(out) :: blparams

        blparams = this%Param

        end subroutine getparam2_Quadrupole

        subroutine getparam3_Quadrupole(this,blength,bnseg,bmapstp,&
                                       btype)
        implicit none
        type (Quadrupole), intent(in) :: this
        double precision, intent(out) :: blength
        integer, intent(out) :: bnseg,bmapstp,btype

        blength = this%Length
        bnseg = this%Nseg
        bmapstp = this%Mapstp
        btype = this%Itype

        end subroutine getparam3_Quadrupole

        subroutine maplinear_Quadrupole(t,tau,xm,this,refpt,Bchg,Bmass)
        implicit none
        include 'mpif.h'
        double precision, intent(in) :: t,tau,Bchg,Bmass
        double precision, dimension(6,6), intent(out) :: xm
        double precision, dimension(6), intent(inout) :: refpt
        type (Quadrupole), intent(in) :: this
        double precision, dimension(18) :: y
        double precision :: h,t0,gi,betai,ui,squi,sqi3
        double precision :: dlti,thli,gf,betaf,uf,squf,sqf3
        double precision :: dltf,thlf,zedge,escale,tfin
        integer  :: mpstp

        xm = 0.0

        zedge = this%Param(1)
        escale = this%Param(2)
        mpstp = this%Mapstp

        y(1)=0.0
        y(2)=0.0
        y(3)=0.0
        y(4)=0.0
        y(5)=refpt(5)
        y(6)=refpt(6)
        y(7)=1.0
        y(8)=0.0
        y(9)=0.0
        y(10)=1.0
        y(11)=1.0
        y(12)=0.0
        y(13)=0.0
        y(14)=1.0
        y(15)=1.0
        y(16)=0.0
        y(17)=0.0
        y(18)=1.0
        h=tau/mpstp
        t0=t
        gi=-refpt(6)
        betai=sqrt(gi**2-1.)/gi
        ui=gi*betai
        squi=sqrt(ui)
        sqi3=sqrt(ui)**3

        dlti=0.0
        thli=0.0

        call rk6i_Quadrupole(h,mpstp,t0,y,18,this,Bchg,Bmass)

        refpt(5)=y(5)
        refpt(6)=y(6)
        gf=-refpt(6)
        betaf=sqrt(gf**2-1.)/gf
        uf=gf*betaf
        squf=sqrt(uf)
        sqf3=sqrt(uf)**3

        tfin = t + tau

        dltf = 0.0
        thlf = 0.0

        xm(1,1)= y(7)*squi/squf+y(9)*squi/squf*dlti
        xm(2,1)=(y(8)-y(7)*dltf)*squi*squf+(y(10)-y(9)*dltf)*squi*squf*&
                 dlti
        xm(1,2)= y(9)/(squi*squf)
        xm(2,2)=(y(10)-y(9)*dltf)*squf/squi
        xm(3,3)= y(11)*squi/squf+y(13)*squi/squf*dlti
        xm(4,3)= &
        (y(12)-y(11)*dltf)*squi*squf+(y(14)-y(13)*dltf)*squi*squf*dlti
        xm(3,4)= y(13)/(squi*squf)
        xm(4,4)=(y(14)-y(13)*dltf)*squf/squi
        xm(5,5)= y(15)*sqi3/sqf3+y(17)*sqi3/sqf3*thli
        xm(6,5)= &
        (y(16)-y(15)*thlf)*sqi3*sqf3+(y(18)-y(17)*thlf)*sqi3*sqf3*thli
        xm(5,6)= y(17)/(sqi3*sqf3)
        xm(6,6)=(y(18)-y(17)*thlf)*sqf3/sqi3
        !print*,"xm: ",xm(1,1),xm(2,1),xm(1,2),xm(2,2),xm(3,3),xm(4,3),xm(4,4)

        end subroutine maplinear_Quadrupole

        subroutine rk6i_Quadrupole(h,ns,t,y,nvar,this,Bchg,Bmass)
        implicit none
        include 'mpif.h'
        integer, intent(in) :: ns,nvar
        double precision, intent(inout) :: t
        double precision, intent(in) :: h,Bchg,Bmass
        double precision, dimension(nvar), intent(inout) :: y
        type (Quadrupole), intent(in) :: this
        double precision, dimension(nvar) :: yt,a,b,c,d,e,f,g,o,p
        double precision:: tint,tt
        integer :: i
        tint=t
        do i=1,ns
          call intfunc1_Quadrupole(t,y,f,this,Bchg,Bmass)
          a=h*f
          yt=y+a/9.d0
          tt=t+h/9.d0
          call intfunc1_Quadrupole(tt,yt,f,this,Bchg,Bmass)
          b=h*f
          yt=y + (a + 3.d0*b)/24.d0
          tt=t+h/6.d0
          call intfunc1_Quadrupole(tt,yt,f,this,Bchg,Bmass)
          c=h*f
          yt=y+(a-3.d0*b+4.d0*c)/6.d0
          tt=t+h/3.d0
          call intfunc1_Quadrupole(tt,yt,f,this,Bchg,Bmass)
          d=h*f
          yt=y + (278.d0*a - 945.d0*b + 840.d0*c + 99.d0*d)/544.d0
          tt=t+.5d0*h
          call intfunc1_Quadrupole(tt,yt,f,this,Bchg,Bmass)
          e=h*f
          yt=y + (-106.d0*a+273.d0*b-104.d0*c-107.d0*d+48.d0*e)/6.d0
          tt = t+2.d0*h/3.d0
          call intfunc1_Quadrupole(tt,yt,f,this,Bchg,Bmass)
          g=h*f
          yt = y+(110974.d0*a-236799.d0*b+68376.d0*c+ &
                  103803.d0*d-10240.d0*e + 1926.d0*g)/45648.d0
          tt = t + 5.d0*h/6.d0
          call intfunc1_Quadrupole(tt,yt,f,this,Bchg,Bmass)
          o=h*f
          yt = y+(-101195.d0*a+222534.d0*b-71988.d0*c-  &
                  26109.d0*d-2.d4*e-72.d0*g+22824.d0*o)/25994.d0
          tt = t + h
          call intfunc1_Quadrupole(tt,yt,f,this,Bchg,Bmass)
          p=h*f
          y = y+(41.d0*a+216.d0*c+27.d0*d+ &
                 272.d0*e+27.d0*g+216.d0*o+41.d0*p)/840.d0
          t=tint+i*h
        enddo

        end subroutine rk6i_Quadrupole

        subroutine intfunc1_Quadrupole(t,y,f,this,Bchg,Bmass)
        implicit none
        include 'mpif.h'
        double precision, intent(in) :: t,Bchg,Bmass
        double precision, dimension(:), intent(in) :: y
        double precision, dimension(:), intent(out) :: f
        type (Quadrupole), intent(in) :: this
        double precision :: gamma0,beta0,gbet,s11,s33,s55
        double precision :: zedge,qmcc,brho,bgrad,zz

        zedge = this%Param(1)
        if(this%Param(3).gt.1.0e-5) then
          zz = t - zedge
          call getfldfrg_Quadrupole(zz,this,bgrad)
        else
          bgrad = this%Param(2)
        endif
        qmcc = Bchg/Bmass

        f(1) = 0.0
        f(2) = 0.0
        f(3) = 0.0
        f(4) = 0.0

        ! synchronous particle:
        gamma0=-y(6)
        beta0=sqrt(gamma0**2-1.)/gamma0
        gbet=sqrt((gamma0-1.0)*(gamma0+1.0))

        f(5)=1.0/(beta0*Scxl)
        f(6)=0.0

        ! matrix elements
        brho=gbet/Clight/qmcc
        s11=bgrad/brho*Scxl
        s33=-bgrad/brho*Scxl
        s55=  0.0

        f(7)=y(8)/Scxl
        f(8)=-s11*y(7)
        f(9)=y(10)/Scxl
        f(10)=-s11*y(9)
        f(11)=y(12)/Scxl
        f(12)=-s33*y(11)
        f(13)=y(14)/Scxl
        f(14)=-s33*y(13)
        f(15)=y(16)/Scxl
        f(16)=-s55*y(15)
        f(17)=y(18)/Scxl
        f(18)=-s55*y(17)

        end subroutine intfunc1_Quadrupole

        !> get external field with displacement and rotation errors.
        subroutine  getflderr_Quadrupole(pos,extfld,this,dx,dy,anglex,&
                                         angley,anglez)
        implicit none
        include 'mpif.h'
        double precision, dimension(4), intent(in) :: pos
        type (Quadrupole), intent(in) :: this
        double precision, intent(in) :: dx,dy,anglex,angley,anglez
        double precision, dimension(6), intent(out) :: extfld
        double precision:: zz,bgrad,zedge
        double precision, dimension(3) :: temp,tmp
        real*8 :: bgradp,bgradpp

        zedge = this%Param(1)
        if(this%Param(3).gt.1.0e-5) then
          zz = pos(3)-zedge
          !call getfldfrg_Quadrupole(zz,this,bgrad)
          !call getfldfrgAna_Quadrupole(zz,this,bgrad)
          call getfldfrgAna2_Quadrupole(zz,this,bgrad,bgradp,bgradpp)
        else
          bgrad = this%Param(2)
        endif

!        dx = this%Param(5)
!        dy = this%Param(6)
!        anglex = this%Param(7)
!        angley = this%Param(8)
!        anglez = this%Param(9)

        temp(1) = pos(1) - dx
        temp(2) = pos(2) - dy
        tmp(1) = temp(1)*cos(anglez) + temp(2)*sin(anglez)
        tmp(2) = -temp(1)*sin(anglez) + temp(2)*cos(anglez)
        tmp(3) = pos(3) - zedge
        temp(1) = tmp(1)*cos(angley)+tmp(3)*sin(angley)
        temp(2) = tmp(2)
        temp(3) = -tmp(1)*sin(angley)+tmp(3)*cos(angley)
        tmp(1) = temp(1)
        tmp(2) = temp(2)*cos(anglex)+temp(3)*sin(anglex)
        tmp(3) = -temp(2)*sin(anglex)+temp(3)*cos(anglex)
        zz = tmp(3)

        extfld(1) = 0.0
        extfld(2) = 0.0
        extfld(3) = 0.0
        extfld(4) = bgrad*tmp(2)
        extfld(5) = bgrad*tmp(1)
        extfld(6) = 0.0

        tmp(1) = extfld(4)
        tmp(2) = extfld(5)*cos(anglex)-extfld(6)*sin(anglex)
        tmp(3) = extfld(5)*sin(anglex)+extfld(6)*cos(anglex)
        temp(1) = tmp(1)*cos(angley)-tmp(3)*sin(angley)
        temp(2) = tmp(2)
        temp(3) = tmp(1)*sin(angley)+tmp(3)*cos(angley)
        extfld(4) = temp(1)*cos(anglez) - temp(2)*sin(anglez)
        extfld(5) = temp(1)*sin(anglez) + temp(2)*cos(anglez)
        extfld(6) = temp(3)

        end subroutine getflderr_Quadrupole

        !> get external field without displacement and rotation errors.
        subroutine getfld_Quadrupole(pos,extfld,this)
        implicit none
        include 'mpif.h'
        double precision, dimension(4), intent(in) :: pos
        type (Quadrupole), intent(in) :: this
        double precision, dimension(6), intent(out) :: extfld
        double precision:: zz,bgrad,zedge,bgradp,bgradpp

        zedge = this%Param(1)
        zz=pos(3)-zedge
        if(this%Param(3).gt.1.0e-5) then
          !call getfldfrg_Quadrupole(zz,this,bgrad)
          !call getfldfrgAna_Quadrupole(zz,this,bgrad)
          call getfldfrgAna2_Quadrupole(zz,this,bgrad,bgradp,bgradpp)
          !bgradp = 0.0
          !bgradpp = 0.0
        else
          bgrad = this%Param(2)
        endif

        extfld(1) = 0.0
        extfld(2) = 0.0
        extfld(3) = 0.0
        if(this%Param(3).gt.0.0) then
          extfld(4) = bgrad*pos(2) - bgradpp*pos(2)**3/6
          extfld(5) = bgrad*pos(1) - bgradpp*pos(1)*pos(2)**2/2
          extfld(6) = bgradp*pos(1)*pos(2)
        else
          extfld(4) = bgrad*pos(2)
          extfld(5) = bgrad*pos(1)
          extfld(6) = 0.0
        endif

        end subroutine getfld_Quadrupole

        !> interpolate the field from the SC rf cavity onto bunch location.
        subroutine getfldfrg_Quadrupole(zz,this,bgrad)
        implicit none
        include 'mpif.h'
        type (Quadrupole), intent(in) :: this
        double precision, intent(in) :: zz
        double precision, intent(out) :: bgrad
        double precision:: hstep,slope
        integer :: klo,khi,k
        integer :: my_rank,ierr

        klo=1
        khi=Ndata
1       if(khi-klo.gt.1) then
          k=(khi+klo)/2
          if(zdat(k).gt.zz)then
             khi=k
          else
             klo=k
          endif
          goto 1
        endif
        hstep=zdat(khi)-zdat(klo)
        slope=(edat(khi)-edat(klo))/hstep
        bgrad =edat(klo)+slope*(zz-zdat(klo))

        end subroutine getfldfrg_Quadrupole

        subroutine getfldfrgAna_Quadrupole(zz,this,bgrad)
        implicit none
        include 'mpif.h'
        type (Quadrupole), intent(in) :: this
        double precision, intent(in) :: zz
        double precision, intent(out) :: bgrad
        double precision :: bb,dd,z1,z2

        bb = this%Param(2)
        z1 = 2*(this%Length - this%Param(3))/2
        z2 = this%Length - z1
        dd = 2*this%Param(4)
        if(zz.lt.z1) then
          bgrad = bb/(1.0+exp(-0.00004+4.518219*(-(zz-z1)/dd-1.5)))
        else if(zz.gt.z2) then
          bgrad = bb/(1.0+exp(-0.00004+4.518219*((zz-z2)/dd-1.5)))
        else
          bgrad = bb
        endif

        !write(1,*)zz,bgrad

        end subroutine getfldfrgAna_Quadrupole

        subroutine getfldfrgAna2_Quadrupole(zz,this,bgrad,bgradp,bgradpp)
        implicit none
        include 'mpif.h'
        type (Quadrupole), intent(in) :: this
        double precision, intent(in) :: zz
        double precision, intent(out) :: bgrad
        double precision :: bgradp,bgradpp
        double precision :: bb,dd,z10,z20,z3,z4,tmpz
        double precision :: c1,c2,s1,s2,s3,s4
        double precision :: tmp1

        c1 = -0.00004d0
        c2 = 4.518219d0
        bb = this%Param(2)
        dd = 2*this%Param(4)
        !coordinate 0 points for Enge function.
        z10 = (this%Length - this%Param(3))/2 !for entrance
        z20 = this%Length - z10               !for exit
        !fringe field range inside coordinate 0 point.
        !this could be different for different Enge coefficients.
        z3 = z10 + 1.5*dd  !for entrance
        z4 = z20 - 1.5*dd  !for exit
        if(zz.lt.z3) then
          tmpz = -(zz - z10)/dd
          s1 = c1 + c2*tmpz
          bgrad = bb/(1.0+exp(s1))
          s2 = -c2/dd !first derivative
          tmp1 = exp(s1)*s2
          bgradp = -bb*tmp1/(1.0+exp(s1))**2
          s3 = 0.0 !second derivative
          s4 = exp(s1)*(s3+s2*s2)
          bgradpp = bb*(-s4/(1.0+exp(s1))**2+2*tmp1*tmp1/(1.0+exp(s1))**3)
        else if(zz.gt.z4) then
          tmpz = (zz - z20)/dd
          s1 = c1 + c2*tmpz
          bgrad = bb/(1.0+exp(s1))
          s2 = c2/dd
          tmp1 = exp(s1)*s2
          bgradp = -bb*tmp1/(1.0+exp(s1))**2
          s3 = 0.0
          s4 = exp(s1)*(s3 + s2*s2)
          bgradpp = bb*(-s4/(1.0+exp(s1))**2+2*tmp1*tmp1/(1.0+exp(s1))**3)
        else
          bgrad = bb
          bgradp = 0.0
          bgradpp = 0.0
        endif

        !write(1,*)zz,bgrad,bgradp,bgradpp

        end subroutine getfldfrgAna2_Quadrupole

        !> quad strength using K
        subroutine transfmapK_Quadrupole(tt,tau,this,refpt,Nplc,pts,qmass)
        implicit none
        include 'mpif.h'
        integer, intent(in) :: Nplc
        double precision, intent(inout) :: tt
        double precision, intent(in) :: tau,qmass
        double precision, dimension(6), intent(inout) :: refpt
        double precision, dimension(6) :: tmp
        type (Quadrupole), intent(in) :: this
        double precision, pointer, dimension(:,:) :: pts
        real*8 :: xm11,xm12,xm21,xm22,xm33,xm34,xm43,xm44,gam,gambetz,&
                  betaz,beta0,rtkstrzz,rtkstr,kstr,gambet0
        integer  :: i
        real*8 :: t,cs,ss

        gambet0 = sqrt(refpt(6)**2-1.0d0)
        beta0 = sqrt(1.0d0-1.0d0/(refpt(6)**2))
        do i = 1, Nplc
          gam = -refpt(6) - pts(6,i)
          gambetz = sqrt(gam**2-1.0d0-pts(2,i)**2-pts(4,i)**2)
          betaz = gambetz/gam
          !each reference particle momentum
          gambetz = sqrt(gam**2-1.0d0)
!          kstr = pts(7,i)*2.997928e8/gambetz*this%Param(2)
!Param(2) is the K defined in MAD, i.e. G/Brho
          kstr = pts(7,i)/qmass*gambet0/gambetz*this%Param(2)
          rtkstr = sqrt(abs(kstr))
          rtkstrzz = rtkstr*tau
          if(kstr.gt.0.0) then
            cs = cos(rtkstrzz)
            ss = sin(rtkstrzz)
            xm11 = cs
            xm12 = ss/rtkstr
            xm21 = -rtkstr*ss
            xm22 = cs
            !xm11 = cos(rtkstrzz)
            !xm12 = sin(rtkstrzz)/rtkstr
            !xm21 = -rtkstr*sin(rtkstrzz)
            !xm22 = cos(rtkstrzz)
            t = exp(rtkstrzz)
            cs= (t*t +1.0d0) / (t + t)
            ss= (t*t -1.0d0) / (t + t)
            xm33 = cs
            xm34 = ss/rtkstr
            xm43 = ss*rtkstr
            xm44 = cs
            !xm33 = cosh(rtkstrzz)
            !xm34 = sinh(rtkstrzz)/rtkstr
            !xm43 = sinh(rtkstrzz)*rtkstr
            !xm44 = cosh(rtkstrzz)
          else if(kstr.lt.0.0) then
            t = exp(rtkstrzz)
            cs= (t*t +1.0d0) / (t + t)
            ss= (t*t -1.0d0) / (t + t)
            xm11 = cs
            xm12 = ss/rtkstr
            xm21 = rtkstr*ss
            xm22 = cs
            !xm11 = cosh(rtkstrzz)
            !xm12 = sinh(rtkstrzz)/rtkstr
            !xm21 = rtkstr*sinh(rtkstrzz)
            !xm22 = cosh(rtkstrzz)
            cs = cos(rtkstrzz)
            ss = sin(rtkstrzz)
            xm33 = cs
            xm34 = ss/rtkstr
            xm43 = -ss*rtkstr
            xm44 = cs
            !xm33 = cos(rtkstrzz)
            !xm34 = sin(rtkstrzz)/rtkstr
            !xm43 = -sin(rtkstrzz)*rtkstr
            !xm44 = cos(rtkstrzz)
          else
            xm11 = 1.0d0
            xm12 = tau
            xm21 = 0.0d0
            xm22 = 1.0d0
            xm33 = 1.0d0
            xm34 = tau
            xm43 = 0.0
            xm44 = 1.0d0
          endif
          tmp(1) = xm11*pts(1,i)+xm12*pts(2,i)/gambetz/Scxl
          tmp(2) = gambetz*Scxl*xm21*pts(1,i)+xm22*pts(2,i)
          tmp(3) = xm33*pts(3,i)+xm34*pts(4,i)/gambetz/Scxl
          tmp(4) = gambetz*Scxl*xm43*pts(3,i)+xm44*pts(4,i)
          tmp(5) = pts(5,i) + (1.0/betaz-1.0/beta0)*tau/Scxl
          tmp(6) = pts(6,i)
          pts(1,i) = tmp(1)
          pts(2,i) = tmp(2)
          pts(3,i) = tmp(3)
          pts(4,i) = tmp(4)
          pts(5,i) = tmp(5)
          pts(6,i) = tmp(6)
        enddo
        refpt(5) = refpt(5) + tau/beta0/Scxl
!        tt = tt + tau

        end subroutine transfmapK_Quadrupole

      end module Quadrupoleclass

max_stars_repo_name: elainejut/IMPACT-Z

max_issues_repo_name: elainejut/IMPACT-Z

max_forks_repo_name: elainejut/IMPACT-Z


text:
      logical function filterWbbmas()
c--- this routine is specific to the "Wbbmas" processes 401-408;
c--F also for processes 520-529
c--- it inspects the jets to check whether an event should be included
c--- routine returns FALSE if event does not pass the process-specific cuts
      implicit none
      include 'constants.f'
      include 'clustering.f'
      include 'jetlabel.f'
      include 'nproc.f'
      logical veto3jets

      if (jets.lt.1) then
         filterWbbmas=.false.
         return
      endif

      filterWbbmas=.true.

c--F  3-jet veto to match with ATLAS analysis. In general we do not want such a veto
      veto3jets = .false.

c---  20 '  f(p1)+f(p2) --> W^+(-->nu(p3)+e^+(p4)) +b(p5)+b~(p6) [massive]'
      if ((nproc .eq. 20) .or. (nproc .eq. 25)) then
        if (jets .lt. 2) then
          filterWbbmas=.false.
          return
c--- note: jets should already contain bq and ba because of bbproc
        endif
      endif

c--  New final state slicing: Wb, W(bb~), Wbb

c--  401 '  f(p1)+f(p2) --> W^+(-->nu(p3)+e^+(p4))+b(p5) [1, 2 or 3 jets, 4FNS]'
      if ((nproc .eq. 401) .or. (nproc .eq. 406)) then
         if ((inclusive .eqv. .false.) .and. (jets .ne. 1)) then
            filterWbbmas=.false.
            return
         endif
         if (jets .eq. 1) then
            if((jetlabel(1).ne.'bq').and.(jetlabel(1).ne.'ba')) then
               filterWbbmas=.false.
               return
            endif
         endif
         if (jets .eq. 2) then
            if ((jetlabel(1).eq.'bb').or.(jetlabel(2).eq.'bb')) then
               filterWbbmas=.false.
               return
            endif
         endif
         if ((jets .eq. 3) .and. (veto3jets .eqv. .true.)) then
            filterWbbmas=.false.
            return
         endif
      endif

c--  402 '  f(p1)+f(p2) --> W^+(-->nu(p3)+e^+(p4))+(b+b~)(p5) [1 or 2 jets, 4FNS]'

      if ((nproc .eq. 402) .or. (nproc .eq. 407)) then
         if ((inclusive .eqv. .false.) .and. (jets .ne. 1)) then
            filterWbbmas=.false.
            return
         endif
         if ((jets .eq. 1) .and. (jetlabel(1) .ne. 'bb')) then
            filterWbbmas=.false.
            return
         endif
         if (jets .eq. 2) then
            if ((jetlabel(1).ne.'bb').and.(jetlabel(2).ne.'bb')) then
               filterWbbmas=.false.
               return
            endif
         endif
         if (jets .eq. 3) then
            filterWbbmas=.false.
            return
         endif
      endif

c--  403 '  f(p1)+f(p2) --> W^+(-->nu(p3)+e^+(p4))+b(p5)+b~(p6) [2 or 3 jets, 4FNS]'
      if ((nproc .eq. 403) .or. (nproc .eq. 408)) then
         if ((jets .eq. 3) .and.
     &       ((inclusive.eqv..false.).or.(veto3jets.eqv..true.))) then
            filterWbbmas=.false.
            return
         endif
         if (jets .eq. 1) then
            filterWbbmas=.false.
            return
         endif
         if (jets .eq. 2) then
            if ((jetlabel(1).eq.'pp').or.(jetlabel(2).eq.'pp')) then
               filterWbbmas=.false.
               return
            endif
         endif
      endif

      return
      end

c--- This code is now deprecated as of v6.1.
c--- The functionality is (mostly) replicated in processes 520-527.

cc--- 420 '  f(p1)+f(p2) --> W^+(-->nu(p3)+e^+(p4))+b(p5) [1 or 2 jets only]'
c      if ((nproc .eq. 420) .or. (nproc .eq. 425)) then
c        if ((jets .gt. 2) .or. (jets .lt. 1)) then
c          filterWbbmas=.false.
c          return
c        endif
c        if ((inclusive .eqv. .false.) .and. (jets .eq. 2)) then
c          filterWbbmas=.false.
c          return
c        endif
c        if (jets .eq. 1) then
c         if ((jetlabel(1) .ne. 'bq') .and. (jetlabel(1) .ne. 'ba')) then
c            filterWbbmas=.false.
c            return
c          endif
c        endif
c        if (jets .eq. 2) then
c          if ((jetlabel(1) .eq. 'bb') .or. (jetlabel(2) .eq. 'bb')) then
c            filterWbbmas=.false.
c            return
c          endif
c        endif
c      endif
c
cc--- 421 '  f(p1)+f(p2) --> W^+(-->nu(p3)+e^+(p4))+(b+b~)(p5) [1 or 2 jets only]'
c      if ((nproc .eq. 421) .or. (nproc .eq. 426)) then
c        if ((jets .gt. 2) .or. (jets .lt. 1)) then
c          filterWbbmas=.false.
c          return
c        endif
c        if ((inclusive .eqv. .false.) .and. (jets .eq. 2)) then
c          filterWbbmas=.false.
c          return
c        endif
c        if (jets .eq. 1) then
c          if ((jetlabel(1) .ne. 'bb')) then
c            filterWbbmas=.false.
c            return
c          endif
c        endif
c        if (jets .eq. 2) then
c         if ((jetlabel(1) .ne. 'bb') .and. (jetlabel(2) .ne. 'bb')) then
c            filterWbbmas=.false.
c            return
c         endif
c        endif
c      endif
c
cc--- 422 '  f(p1)+f(p2) --> W^+(-->nu(p3)+e^+(p4))+b(p5)+jet(p6) [2 or 3 jets only]'
c      if ((nproc .eq. 422) .or. (nproc .eq. 427)) then
c        if (jets .lt. 2) then
c          filterWbbmas=.false.
c          return
c        endif
c        if ((inclusive .eqv. .false.) .and. (jets .eq. 3)) then
c          filterWbbmas=.false.
c          return
c        endif
c        if (jets .eq. 2) then
c          if ((jetlabel(1) .eq. 'bb') .or. (jetlabel(2) .eq. 'bb')) then
c            filterWbbmas=.false.
c            return
c          endif
c        endif
c      endif

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


MODULE cri_mosaic_8bin_aq_Model

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  Completely defines the model cri_mosaic_8bin_aq
!    by using all the associated modules
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  USE cri_mosaic_8bin_aq_Precision
  USE cri_mosaic_8bin_aq_Parameters
  USE cri_mosaic_8bin_aq_Global
  USE cri_mosaic_8bin_aq_Function
  USE cri_mosaic_8bin_aq_Integrator
  USE cri_mosaic_8bin_aq_Rates
  USE cri_mosaic_8bin_aq_Jacobian
  USE cri_mosaic_8bin_aq_LinearAlgebra
  USE cri_mosaic_8bin_aq_Monitor
  USE cri_mosaic_8bin_aq_Util

END MODULE cri_mosaic_8bin_aq_Model


text:
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!                                                                      %
! Copyright (C) 2015 Adrian Martinez Vargas                            %
!                                                                      %
! This software may be modified and distributed under the terms        %
! of the MIT license.  See the LICENSE.txt file for details.           %
!                                                                      %
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!-----------------------------------------------------------------------
!
! The subroutine below isbased on GSLIB code,
! version 2.0 of the gslib code written in fortran 77
!
!-----------------------------------------------------------------------

!*********************************************************************************
!     Subroutines in GSLIB (auxiliary functions)
!*********************************************************************************
include 'setrot.f90'
include 'cova3.f90'
include 'sqdist.f90'
! include 'getindx.f90'
! include 'locate.f90'
! this includes sortem and dsortem
! include 'sortem.f90'
! include 'acorni.f90'
! include 'gauinv.f90'
! include 'powint.f90'
! include 'gcum.f90'

subroutine block_covariance(xdb,ydb,zdb, ndb, &
                            nst,it,c0,cc,aa,aa1,aa2,ang1,ang2,ang3, &
                            unbias,cbb)
    !-----------------------------------------------------------------------

    !                    Block Covariance
    !                    *****************************

    ! This subroutine calculate the block covariance associated with a variogram
    ! model specified by a nugget effect and nested varigoram structures.
    ! The anisotropy definition can be different for each nested structure.
    ! The block size is defined with input discretization points with
    ! arbitrary locations (for example regular discretization, random
    ! discretization or discretization points in an irregular polygon)

    ! INPUT VARIABLES:

    !   xdb,ydb,zdb      coordinates of discretization points
    !   ndb              number of discretization points
    !   nst(ivarg)       number of nested structures (maximum of 4)
    !   c0(ivarg)        isotropic nugget constant
    !   it(i)            type of each nested structure:
    !                      1. spherical model of range a;
    !                      2. exponential model of parameter a;
    !                           i.e. practical range is 3a
    !                      3. gaussian model of parameter a;
    !                           i.e. practical range is a*sqrt(3)
    !                      4. power model of power a (a must be gt. 0  and
    !                           lt. 2).  if linear model, a=1,c=slope.
    !                      5. hole effect model
    !   cc(i)            multiplicative factor of each nested structure.
    !                      (sill-c0) for spherical, exponential,and gaussian
    !                      slope for linear model.
    !   aa,aa1,aa2       parameters "a" of each nested structure.
    !   aa,aa1,aa2       rotation angles for each nested structure
    !
    ! OUTPUT VARIABLES:

    !   unbias           unbias variable, need internally for kriging
    !   cbb              block covariance

    ! Notes

    ! This functions was created from code lines in kt3d program
    !
    !
    ! Adrian Martinez  2015
    !-----------------------------------------------------------------------

    IMPLICIT NONE

    ! in
    ! variogram
    integer, intent(in)                 :: nst
    real*8, intent(in),   dimension(1)  :: c0
    integer, intent(in), dimension(nst) :: it
    real*8, intent(in), dimension(nst)  :: cc,aa, aa1, aa2,ang1,ang2,ang3
    !discretization points of the last block
    integer, intent (in) :: ndb   ! number of discretization points
    real*8, intent (in), dimension(ndb) :: xdb,ydb,zdb      !coordinates of discretization points

    !out
    real*8, intent(out) :: unbias, cbb

    ! internal
    ! variogram
    real*8 :: cmax,covmax, cov
    real*8, dimension(nst) :: anis1, anis2
    real*8, dimension(nst,3,3) :: rotmat
    real*8 ::  EPSLON=1.e-20, PMX=999.
    integer :: i, is, j

    do i=1,nst
        anis1(i) = aa1(i) / max(aa(i),EPSLON)
        anis2(i) = aa2(i) / max(aa(i),EPSLON)

        ! print *, aa(i), aa1(i), aa2(i), anis1(i), anis2(i)

        if(it(i).eq.4) then
              if(aa(i).lt.0.0) stop ' INVALID power variogram'
              if(aa(i).gt.2.0) stop ' INVALID power variogram'
        end if
    end do

    ! get the rotation matrix
    do is=1,nst
        ! bug fixed: setrot with doubles require dsetrot
        call dsetrot(ang1(is),ang2(is),ang3(is),anis1(is),anis2(is), &
        is,nst,rotmat)
        if(it(is) == 4) then
            covmax = covmax + PMX
        else
            covmax = covmax + cc(is)
        endif
    end do

    ! Calculate Block Covariance. Check for point kriging.
    ! bug fixed: dcova3 with doubles require dsetrot
    call dcova3(xdb(1),ydb(1),zdb(1),xdb(1),ydb(1),zdb(1),1,nst, &
                c0,it,cc,aa,1,nst,rotmat,cmax,cov)

    ! Set the 'unbias' variable so that the matrix solution is more stable

    unbias = cov
    cbb    = dble(cov)
    if(ndb > 1) then
        cbb = 0.0
        do i=1,ndb
            do j=1,ndb
                call dcova3(xdb(i),ydb(i),zdb(i),xdb(j),ydb(j),zdb(j), &
                           1,nst,c0,it,cc,aa,1,nst,rotmat,cmax,cov)
                if(i == j) cov = cov - c0 (1)
                cbb = cbb + dble(cov)
            end do
        end do
        cbb = cbb/dble(ndb*ndb)
    end if

    return

end subroutine block_covariance

max_stars_repo_name: cmrajan/pygslib

max_issues_repo_name: kaufmanno/pygslib

max_forks_repo_name: kaufmanno/pygslib


text:
PROGRAM test_abstraction2
 USE mo_column , ONLY: compute_column
 REAL :: q ( 1 : 10 , 1 : 10 , 1 : 60 )
 REAL :: t ( 1 : 10 , 1 : 10 , 1 : 60 )
 INTEGER :: nx
 INTEGER :: ny
 INTEGER :: nz
 INTEGER :: i
 INTEGER :: j

 nx = 10
 ny = 10
 nz = 60
 DO i = 1 , nx , 1
  DO j = 1 , ny , 1
   q ( i , j , 1 ) = 0.0
   t ( i , j , 1 ) = 0.0
  END DO
 END DO
!$acc data pcreate(q(:,:,:),t(:,:,:))
!$acc update device(q(:,:,:),t(:,:,:))
 CALL compute_column ( nz , q ( : , : , : ) , t ( : , : , : ) , nx = nx , ny =&
  ny )
!$acc update host(q(:,:,:),t(:,:,:))
!$acc end data
 PRINT * , sum ( q )
 PRINT * , sum ( t )
END PROGRAM test_abstraction2

max_stars_repo_name: peclatj/claw-compiler

max_issues_repo_name: peclatj/claw-compiler

max_forks_repo_name: peclatj/claw-compiler


text: Davis relevant clipart and general use photos for illustration of entries and Events Board events.  These are (like everything on the wiki) all CCBy, so feel free to use them with attribution on fliers and such.  Ideally, these are modified for conceptual use... images specific to a topic should be used directly on the relevant entry.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE POLARC(L,CBUF,AMP,PHZ)
C     POLARC COMPUTES THE POLAR COORDINATES FROM RECTANGULAR COORDINATES.
C
C ARGUMENTS:
C L    - THE NUMBER OF COORDINATES TO CONVERT.
C        THE NUMBER OF COMPLEX PAIRS IN ARRAY CBUF.
C CBUF  - AN ARRAY OF COMPLEX RECTANGULAR COORDINATES.
C AMP  - THE ARRAY OF AMPLITUDE COMPONENTS CALCULATED BY POLAR.
C PHZ  - THE ARRAY OF PHASES CALCULATED.
C
C   TAKEN FROM ENDERS ROBINSON'S BOOK "MULTICHANNEL TIME SERIES ANALYSIS
C  WITH DIGITAL COMPUTER PROGRAMS": HOLDEN-DAY, 1967, SAN FRANSISCO.
C
c  mod 1 Apr 10 - remove arithmetic IF statements
c
      DIMENSION CBUF(1),AMP(L),PHZ(L)
      COMPLEX CBUF
      PI=3.14159265
      DO 110 I=1,L
      RE=REAL(CBUF(I))
      XIM=AIMAG(CBUF(I))
      AMP(I)=SQRT(RE*RE+XIM*XIM)
c      IF(XIM) 10,20,30
      IF( xim .EQ. 0 ) GOTO 20
      IF( xim .GT. 0 ) GOTO 30
c10    IF(RE) 40,50,60
      IF( re .LT. 0 ) GOTO 40
      IF( re .EQ. 0 ) GOTO 50
      IF( re .GT. 0 ) GOTO 60
c20    IF(RE) 70,80,60
20    IF( re .LT. 0 ) GOTO 70
      IF( re .EQ. 0 ) GOTO 80
      IF( re .GT. 0 ) GOTO 60
c30    IF(RE) 90,100,60
30    IF( re .LT. 0 ) GOTO 90
      IF( re .EQ. 0 ) GOTO 100
      IF( re .GT. 0 ) GOTO 60
40    PHZ(I)=ATAN(XIM/RE)-PI
      GO TO 110
50    PHZ(I)=-PI/2.
      GO TO 110
60    PHZ(I)=ATAN(XIM/RE)
      GO TO 110
70    PHZ(I)=-PI
      GO TO 110
80    PHZ(I)=0.
      GO TO 110
90    PHZ(I)=ATAN(XIM/RE)+PI
      GO TO 110
100   PHZ(I)=PI/2.
110   CONTINUE
      RETURN
      END

max_stars_repo_name: apthorpe/sioseis

max_issues_repo_name: henkart/sioseis

max_forks_repo_name: henkart/sioseis


SUBROUTINE gser(gamser,a,x,gln)
INTEGER*4 ITMAX
REAL*8 a,gamser,gln,x,EPS
PARAMETER (ITMAX=100,EPS=3.d-7)
!CU    USES gammln
INTEGER*4 n
REAL*8 ap,del,sum,gammln
EXTERNAL gammln
gln=gammln(a)
if(x.le.0.)then
	if(x.lt.0.)pause 'x < 0 in gser'
	gamser=0.
	return
endif
ap=a
sum=1./a
del=sum
do n=1,ITMAX
	ap=ap+1.
	del=del*x/ap
	sum=sum+del
	if(abs(del).lt.abs(sum)*EPS)goto 1
end do
pause 'a too large, ITMAX too small in gser'
1     gamser=sum*exp(-x+a*log(x)-gln)
return
END


text:
PROGRAM F014

    ! Copyright 2021 Melwyn Francis Carlo

    IMPLICIT NONE

    INTEGER :: N_COUNT = 1E5

    INTEGER :: MAX_COUNT        = 0
    INTEGER :: MAX_COUNT_NUMBER = 0

    INTEGER :: I_COUNT

    INTEGER (KIND=8) :: N

    DO WHILE (N_COUNT < 1E6)

        I_COUNT = 0
        N = N_COUNT;

        DO WHILE (N /= 1)

            IF (MOD(N, INT(2, 8)) == 0) THEN
                N = N / 2;
            ELSE
                N = (3 * N) + 1;
            END IF

            I_COUNT = I_COUNT + 1

        END DO

        I_COUNT = I_COUNT + 1

        IF (I_COUNT > MAX_COUNT) THEN
            MAX_COUNT        = I_COUNT;
            MAX_COUNT_NUMBER = N_COUNT;
        END IF

        N_COUNT = N_COUNT + 1

    END DO

    PRINT ('(I0)'), MAX_COUNT_NUMBER

END PROGRAM F014

max_stars_repo_name: melwyncarlo/ProjectEuler

max_issues_repo_name: melwyncarlo/ProjectEuler

max_forks_repo_name: melwyncarlo/ProjectEuler


text:
! { dg-do run }
! PR24268 Test case derived from example given by Iwan Kawrakow
! Embedded spaces in format strings should be ignored.
! Prepared by Jerry DeLisle  <jvdelisle@gcc.gnu.org>
      program pr24268
      real x
      character*13 line
      line = "12.34"
      read(line,*) x
      write(line,10) x
 10                                                            format(g1
     *      1.4)
      if (line.ne."  12.34") call abort()
      line = ""
      write(line,20) x
 20   format(t r 2 , g 1 1 . 4)
      if (line.ne."    12.34") call abort()
      end

max_stars_repo_name: vidkidz/crossbridge

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


text:
****************************************************************

      SUBROUTINE tred2(a,n,np,d,e)
      INTEGER n,np
      DOUBLE PRECISION a(np,np),d(np),e(np)
      INTEGER i,j,k,l
      DOUBLE PRECISION f,g,h,hh,scale
      do 18 i=n,2,-1
        l=i-1
        h=0.d0
        scale=0.d0
        if(l.gt.1)then
          do 11 k=1,l
            scale=scale+abs(a(i,k))
11        continue
          if(scale.eq.0.d0)then
            e(i)=a(i,l)
          else
            do 12 k=1,l
              a(i,k)=a(i,k)/scale
              h=h+a(i,k)**2
12          continue
            f=a(i,l)
            g=-sign(sqrt(h),f)
            e(i)=scale*g
            h=h-f*g
            a(i,l)=f-g
            f=0.d0
            do 15 j=1,l
C     Omit following line if finding only eigenvalues
              a(j,i)=a(i,j)/h
              g=0.d0
              do 13 k=1,j
                g=g+a(j,k)*a(i,k)
13            continue
              do 14 k=j+1,l
                g=g+a(k,j)*a(i,k)
14            continue
              e(j)=g/h
              f=f+e(j)*a(i,j)
15          continue
            hh=f/(h+h)
            do 17 j=1,l
              f=a(i,j)
              g=e(j)-hh*f
              e(j)=g
              do 16 k=1,j
                a(j,k)=a(j,k)-f*e(k)-g*a(i,k)
16            continue
17          continue
          endif
        else
          e(i)=a(i,l)
        endif
        d(i)=h
18    continue
C     Omit following line if finding only eigenvalues.
      d(1)=0.d0
      e(1)=0.d0
      do 24 i=1,n
C     Delete lines from here ...
        l=i-1
        if(d(i).ne.0.d0)then
          do 22 j=1,l
            g=0.d0
            do 19 k=1,l
              g=g+a(i,k)*a(k,j)
19          continue
            do 21 k=1,l
              a(k,j)=a(k,j)-g*a(k,i)
21          continue
22        continue
        endif
C     ... to here when finding only eigenvalues.
        d(i)=a(i,i)
C     Also delete lines from here ...
        a(i,i)=1.d0
        do 23 j=1,l
          a(i,j)=0.d0
          a(j,i)=0.d0
23      continue
C     ... to here when finding only eigenvalues.
24    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 5.d0).

C***************************************************************

      SUBROUTINE tqli(d,e,n,np,z)
      INTEGER n,np
      DOUBLE PRECISION d(np),e(np),z(np,np)
CU    USES pythag
      INTEGER i,iter,k,l,m
      DOUBLE PRECISION b,c,dd,f,g,p,r,s,pythag
      do 11 i=2,n
        e(i-1)=e(i)
11    continue
      e(n)=0.d0
      do 15 l=1,n
        iter=0
1       do 12 m=l,n-1
          dd=abs(d(m))+abs(d(m+1))
          if (abs(e(m))+dd.eq.dd) goto 2
12      continue
        m=n
2       if(m.ne.l)then
*         if(iter.eq.30)pause 'too many iterations in tqli'
          iter=iter+1
          g=(d(l+1)-d(l))/(2.d0*e(l))
          r=pythag(g,1.d0)
          g=d(m)-d(l)+e(l)/(g+sign(r,g))
          s=1.d0
          c=1.d0
          p=0.d0
          do 14 i=m-1,l,-1
            f=s*e(i)
            b=c*e(i)
            r=pythag(f,g)
            e(i+1)=r
            if(r.eq.0.d0)then
              d(i+1)=d(i+1)-p
              e(m)=0.d0
              goto 1
            endif
            s=f/r
            c=g/r
            g=d(i+1)-p
            r=(d(i)-g)*s+2.d0*c*b
            p=s*r
            d(i+1)=g+p
            g=c*r-b
C     Omit lines from here ...
            do 13 k=1,n
              f=z(k,i+1)
              z(k,i+1)=s*z(k,i)+c*f
              z(k,i)=c*z(k,i)-s*f
13          continue
C     ... to here when finding only eigenvalues.
14        continue
          d(l)=d(l)-p
          e(l)=g
          e(m)=0.d0
          goto 1
        endif
15    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 5.d0).

C***************************************************************

      FUNCTION pythag(a,b)
      DOUBLE PRECISION a,b,pythag
      DOUBLE PRECISION absa,absb
      absa=abs(a)
      absb=abs(b)
      if(absa.gt.absb)then
        pythag=absa*sqrt(1.d0+(absb/absa)**2)
      else
        if(absb.eq.0.d0)then
          pythag=0.d0
        else
          pythag=absb*sqrt(1.d0+(absa/absb)**2)
        endif
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 5.d0).

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


! chk_parity.f90 --
!     Check if the compiler supports the Fortran 2008 function parity
!
program chk_parity
    implicit none

    write( *, '(a,l)' ) 'The parity of [True,False,True] is ', parity((/.true.,.false.,.true./))
end program chk_parity


text:
	SUBROUTINE GQGGRF  ( ixtyp, iytyp, kx, ky, xll, yll, xur,
     +			     yur,   iret )
C************************************************************************
C* GQGGRF								*
C* 									*
C* This subroutine returns the current coordinate system definition	*
C* for a grid which is evenly spaced in a graph coordinate system.	*
C* The grid coordinate system is defined by GSGGRF.			*
C* 									*
C* GQGGRF  ( IXTYP, IYTYP, KX, KY, XLL, YLL, XUR, YUR, IRET )		*
C*									*
C* Output parameters:							*
C*	IXTYP		INTEGER		X coordinate type 		*
C*					  1 = linear			*
C*					  2 = logarithmic		*
C*					  3 = ** KAPPA (2/7)		*
C*					  5 = polar (R)			*
C*	IYTYP		INTEGER		Y coordinate type		*
C*					  1 = linear			*
C*					  2 = logarithmic		*
C*					  3 = ** KAPPA (2/7)		*
C*					  5 = polar (THETA)		*
C*	KX		INTEGER		Number of x grid points		*
C*	KY		INTEGER		Number of y grid points		*
C*	XLL 		REAL		Lower left X value		*
C*	YLL		REAL		Lower left Y value		*
C*	XUR		REAL		Upper right X value		*
C*	YUR		REAL		Upper right Y value		*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* J. Vilardo/RDS	11/84	GEMPLT Version 3.0			*
C* M. Goodman/RDS	 5/85	GEMPLT Version 3.1			*
C* M. desJardins/GSFC	 6/88	Changed calling sequence to have kx,ky	*
C* L. Williams/EAi       3/94   Removed blank comments from header      *
C* A. Hardy/GSC          6/98   Cleaned up prolog                       *
C************************************************************************
	INCLUDE		'FUNCCODE.PRM'
	INCLUDE		'ERROR.PRM'
C
	INTEGER		isend (4)
	REAL		rrcv (4)
C------------------------------------------------------------------------
C*	Load input parameters into buffer then write them to the mailbox.
C
	isend (1) = 2
	isend (2) = FQGGRF
C
	CALL GPUT  ( isend, 2, iret  )
	IF  ( iret .ne. NORMAL )  RETURN
C
	CALL GGET  ( iret, 1, ierr )
	IF  ( ierr .ne. NORMAL )  THEN
	    iret = ierr
	    RETURN
	END IF
C
	CALL GGET  ( isend, 4, ierr )
	IF  ( ierr .eq. NORMAL )  THEN
	    ixtyp  = isend ( 1 )
	    iytyp  = isend ( 2 )
	    kx     = isend ( 3 )
	    ky     = isend ( 4 )
	  ELSE
	    iret   = ierr
	    RETURN
	END IF
C
	CALL GGETR  ( rrcv, 4, ierr )
	IF  ( ierr .eq. NORMAL )  THEN
	    xll    = rrcv ( 1 )
	    yll    = rrcv ( 2 )
	    xur    = rrcv ( 3 )
	    yur    = rrcv ( 4 )
	  ELSE
	    iret   = ierr
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
C> \file bl.f
C> \brief Solve the Blasius equation using spectral collocation
C> \author S. Scott Collis

c******************************************************************************
      program Solve_Similar_BL
c******************************************************************************
c
c     Purpose:  This program solves the Boundary Layer similarity equation
c               using a Chebyshev collocation method.  Written originally
c               for ME 308 and ME 351B final projects.
c
c               This code works reasonably well but the convergence is rather
c               slow.  I think it can be made faster but I am not sure that
c               it is worth it.  I now have an RK4 boundary layer equation
c               integrator that works quite well.  See CONTEBL.F for example.
c
c     Author:   S. Scott Collis
c
c     Date:     2-22-92
c
c     Revised:  9-18-92
c
c******************************************************************************
      integer     idim
      parameter   (idim=128)
      real        u(0:idim), v(0:idim), xi(0:idim), th(0:idim)
      real        maxerr
      integer     i, j, k, n, nbig, maxcount

      real        gamma, etaout, phi1(0:idim), phi2(0:idim), beta, alpha
      common      /metrics/  gamma, etaout, phi1, phi2, beta, alpha

      write (*,10)
 10   format (/,/,10x,'Solve Boundary Layer Similarity Equation')
      write (*,20)
 20   format (/,1x,'Enter the number of modes ==> ',$)
      read (*,*) n
      write (*,60)
 60   format (/,1x,'Enter Beta ==> ',$)
      read (*,*) beta
      write (*,30)
 30   format (/,1x,'Enter alpha (alpha < 0) ==> ',$)
      read (*,*) alpha
      write (*,40)
 40   format (/,1x,'Enter the max allowable error ==> ',$)
      read (*,*) maxerr
      write (*,70)
 70   format (/,1x,'Enter the max number of iterations ==> ',$)
      read (*,*) maxcount
      write (*,50)
 50   format (/,1x,'Interpolate to mesh size ==> ',$)
      read (*,*) nbig

      gamma = 1.2
      etaout = 15.0

      call MAKE_MESH (n, th, xi)
      call MAKE_METRICS (n, xi)
      call INIT_PROFILE (n, u, v, th, xi)
      call SOLVE_BL (n, u, v, th, xi, maxerr, maxcount)
      call SAVE_VELOCITY(n,u,v,xi)
      call CHEBYINT(n,u,v,nbig)

      stop
      end
C******************************************************************************
      subroutine MAKE_MESH(n, th, xi)
C******************************************************************************
C
C     Setup the collocation points in the mapped coordinate, xi and in
C     chebyshev space, th.
C
C******************************************************************************
      integer    n
      real       th(0:n), xi(0:n)
      real       pi, dth
      integer    i

      pi = Acos (-1.0)
      dth = pi/float(n)

      do i = 0, n
        th(i) = i*dth
        xi(i) = cos(th(i))
      end do

      return
      end
C******************************************************************************
      subroutine MAKE_METRICS(n, xi)
C******************************************************************************
C
C     Compute the first and second metrics for the coordinate transformation
C     I am using the truncated mapping of Street et. al.
C
C******************************************************************************
      integer    n
      real       xi(0:n), A, B
      integer    i

      parameter   (idim=128)
      real        gamma, etaout, phi1(0:idim), phi2(0:idim), beta, alpha
      common      /metrics/  gamma, etaout, phi1, phi2, beta, alpha

      write (*,*)
      write (*,*) 'Mesh and Metrics'
      write (*,*)
      do i = 0, n
        A = etaout/2.*(1-TANH(gamma))
        B = gamma/2.
        eta = etaout*(1-tanh(gamma))*.5*(xi(i)+1)/
     .               (1-tanh(gamma/2.*(xi(i)+1)))
        phi1(i) = 1./(A/2.*(1+EXP(2.*B*(1+xi(i))))+A*B*
     .            EXP(2.*B*(1+xi(i)))*(1+xi(i)))
        phi2(i) = -2.*A*B*(1+B*(xi(i)+1))*EXP(2.*B*(xi(i)+1))*phi1(i)**3
        write (*,10) eta, xi(i), phi1(i), phi2(i)
  10    format (4(es12.5,1x))
      end do

      return
      end
C******************************************************************************
      subroutine INIT_PROFILE(n, u, v, th, xi)
C******************************************************************************
C
C     Setup the initial boundary layer profiles
C
C******************************************************************************
      integer     n
      real        u(0:n), v(0:n), th(0:n), xi(0:n)
      real        eta
      integer     i, ans

      parameter   (idim=128)
      real        gamma, etaout, phi1(0:idim), phi2(0:idim), beta, alpha
      common      /metrics/  gamma, etaout, phi1, phi2, beta, alpha

      write (*,10)
  10  format (/,1x,'Read in a profile (1,0) ==> ',$)
      read (*,*) ans

      if (ans.eq.1) then
        call READ_VELOCITY(n,u,v)
      else
        do i = 0, n
          eta = etaout*(1-tanh(gamma))*.5*(xi(i)+1)/
     .                 (1-tanh(gamma/2.*(xi(i)+1)))
          u(i) = .5*(xi(i)+1)
        end do
c
c       Now given the initial u, v must be set to satisfy continuity
c
        call COMPUTE_V (n, u, v, xi)
      end if

      write (*,*) 'Initial profiles'
      write (*,*)
      do i = 0, n
        eta = etaout*(1-tanh(gamma))*.5*(xi(i)+1)/
     .               (1-tanh(gamma/2.*(xi(i)+1)))
        write (*,12) eta, u(i), v(i)
  12    format (1x,3(e12.5,2x))
      end do
      write (*,*)

      return
      end
C******************************************************************************
      subroutine SPECTRAL(n, u, v, th, xi, L)
C******************************************************************************
C
C     Compute the spectral operator given a velocity field using Fast
C     Chebyshev transform
C
C******************************************************************************
      integer     n
      real        u(0:n), v(0:n), th(0:n), xi(0:n), L(1:n-1)
      integer     i, j
      real        eta

      parameter   (idim=128)
      real        gamma, etaout, phi1(0:idim), phi2(0:idim), beta, alpha
      common      /metrics/  gamma, etaout, phi1, phi2, beta, alpha

      real        u1(0:idim), u2(0:idim)

      if (n.gt.idim) then
        write (*,*) 'Error.  N > idim in SPECTRAL'
        stop
      end if
c
c     Compute the first derivative
c
      call DCHEBYSHEV (u, u1, u2, n)

c      write (*,*)
c      do i = 0, n
c        eta = etaout*(1-tanh(gamma))*.5*(xi(i)+1)/
c     .               (1-tanh(gamma/2.*(xi(i)+1)))
c        write (*,10) eta,u(i),phi1(i)*u1(i),
c     .             (phi1(i)**2*u2(i)+phi2(i)*u1(i)),phi1(i),phi2(i)
c   10   format (1x,6(e11.5,1x))
c      end do
c
C     Now form the spectral operator, L
c
      do i = 1, n-1
        L(i) = phi1(i)**2*u2(i)+(phi2(i)-v(i)*phi1(i))*u1(i)
      end do

      return
      end

C******************************************************************************
      subroutine COMPUTE_V (n, u, v, xi)
C******************************************************************************
C
C     Compute v(i) from continuity given u(i).
C
C******************************************************************************
      integer     n, i
      real        u(0:n), v(0:n), xi(0:n)

      parameter   (idim=128)
      real        gamma, etaout, phi1(0:idim), phi2(0:idim), beta, alpha
      common      /metrics/  gamma, etaout, phi1, phi2, beta, alpha

      real        utemp(0:idim)

      do i = 0, n
        utemp(i) = -1./phi1(i)*u(i)
      end do
      call ICHEBYSHEV (utemp,v,n)

      return
      end

C******************************************************************************
      subroutine FINITE_DIFF(n, u, v, xi, Ha, Hb, Hc, b, L)
C******************************************************************************
C
C     Generate the second order finite difference preconditioning matrix.
C     Note that central difference is used which results in a tridiagonal
C     matrix.  Note that phi1 and phi2 are the first and
C     second metrics of the coordinate transformation.
C
C******************************************************************************
      integer     n
      real        u(0:n), v(0:n), xi(0:n), Ha(n-1), Hb(n-1), Hc(n-1)
      real        b(1:n-1), L(1:n-1)
      real        x10, x11, x01
      integer     i, j

      parameter   (idim=128)
      real        gamma, etaout, phi1(0:idim), phi2(0:idim), beta, alpha
      common      /metrics/  gamma, etaout, phi1, phi2, beta, alpha
c
c     Finite difference operator must also have beta(1-u**2) term. maybe
c
      do i = 1, n-1
        x10 = xi(i+1)-xi(i)
        x11 = xi(i+1)-xi(i-1)
        x01 = xi(i)-xi(i-1)
c
c       Richardson Iteration
c
        Ha(i) = (phi1(i)**2/x01-(phi2(i)-v(i)*phi1(i)))/x11
        Hb(i) = -1.0*phi1(i)**2*(1./x10 + 1./x01)/x11
        Hc(i) = (phi1(i)**2/x10+(phi2(i)-v(i)*phi1(i)))/x11
c
c       Simple pseudo time-step
c
c       Ha(i) = -(phi1(i)**2/x01-(phi2(i)-v(i)*phi1(i)))/x11
c       Hb(i) =  alpha+phi1(i)**2*(1./x10 + 1./x01)/x11
c       Hc(i) = -(phi1(i)**2/x10+(phi2(i)-v(i)*phi1(i)))/x11
      end do
      Ha(1)   = 0.0
      Hc(n-1) = 0.0
c
c     Form the R.H.S. and include the boundary conditions.  Note that
c     u(0) = 1.0 and u(n) = 0.0
c
      do i = 1, n-1
        b(i) = alpha*(L(i)+beta*(1.-u(i)**2))
      end do
c      i = 1
c      x10 = xi(i+1)-xi(i)
c      x11 = xi(i+1)-xi(i-1)
c      x01 = xi(i)-xi(i-1)
c      b(i) = b(i)-(phi1(i)**2/x01-(phi2(i)-v(i)*phi1(i)))/x11*u(0)
c      i = n
c      x10 = xi(i+1)-xi(i)
c      x11 = xi(i+1)-xi(i-1)
c      x01 = xi(i)-xi(i-1)
c      b(i) = b(i)-(phi1(i)**2/x10+(phi2(i)-v(i)*phi1(i)))/x11*u(n)

      return
      end

C******************************************************************************
      subroutine SOLVE_BL(n, u, v, th, xi, maxerr, maxcount)
C******************************************************************************
C
C     Solve the similar boundary layer equations using a Chebyshev collocation
C     method.  On input u, and v contain the initial 'guess' and the non-linear
C     problem is solved iteratively starting from these inital profiles.  On
C     output u and v contain the solution
C
C******************************************************************************
      integer     n
      real        u(0:n), v(0:n), th(0:n), xi(0:n), maxerr
      parameter   (idim=128)
      real        du(idim), b(idim), L(idim)
      real        Ha(idim), Hb(idim), Hc(idim)
      real        residual, eta, oldresidual
      integer     i, j, icount, maxcount

      real        gamma, etaout, phi1(0:idim), phi2(0:idim), beta, alpha
      common      /metrics/  gamma, etaout, phi1, phi2, beta, alpha

      real        H(idim,idim), work(idim), cond, z(idim), oa(idim,idim)
      integer     ipvt(idim)
      real        Ldiff
      real        u1(0:idim), u2(0:idim)

      if (n.gt.idim) then
        write (*,*) 'Error:  N > idim in SOLVE_BL'
        stop
      end if

      write(*,*) "Beginning Iterative Solve"
      write(*,*)

      residual = 100.
      oldresidual = 1000.
      icount = 0
      do while ((residual .gt. maxerr).and.(icount.lt.maxcount))
        icount = icount + 1
c        do i = 1, idim
c          du(i) = 0.0
c          Ha(i) = 0.0
c          Hb(i) = 0.0
c          Hc(i) = 0.0
c          L(i) = 0.0
c          b(i) = 0.0
c        end do
c
c       Setup the spectral R.H.S.
c
        call SPECTRAL (n, u, v, th, xi, L)
c
c       Setup the finite difference opperator
c
        call FINITE_DIFF (n, u, v, xi, Ha, Hb, Hc, b, L)
c
c       Build the matrix
c
c        do i = 2, n-2
c          H(i,i-1) = Ha(i)
c          H(i,i)   = Hb(i)
c          H(i,i+1) = Hc(i)
c        end do
c        H(1,1) = Hb(1)
c        H(1,2) = Hc(1)
c        H(n-1,n-2) = Ha(n-1)
c        H(n-1,n-1) = Hb(n-1)
c        ndim = idim
c
c       H(u) should about equal L
c
c        write (*,*)
c        write (*,*) 'Spectral Finite Difference Comparison'
c        write (*,*)
c        do i = 1, n-1
c          Ldiff = 0.0
c          do j = 1, n-1
c            Ldiff = Ldiff + H(i,j)*u(j)
c          end do
c          if (i.eq.1) then
c            Ldiff = Ldiff - b(1) + L(i)+beta*(1.-u(i)**2)
c          end if
c          write (*,*) xi(i), Ldiff, L(i)
c        end do
c        write (*,*)

c        call DECOMP(ndim,n-1,H,COND,IPVT,WORK,OA,Z)
c        write (*,*) 'Condition = ',cond
c        call SOLVE(ndim,n-1,H,B,IPVT,du)
c
        call TRIDIAG (n-1,du,Ha,Hb,Hc,b)
c
c        try without preconditioning
c
c         do i = 1, n-1
c           du(i) = b(i)
c         end do
c
c       Update u vector, of course the b.c. don't change
c
        do i = 1, n-1
          u(i) = u(i) + du(i)
        end do
c
c       Now Chebyshev integrate to get the new v
c
        call COMPUTE_V (n, u, v, xi)
c
c       Check the residual to seen if it is small enough
c
        oldresidual = residual
        residual = 0.0
        do i = 1, n-1
          residual = residual + (du(i))**2
        end do
        residual = SQRT(residual/(n-1))
        write(*,*) icount, residual
        if ((mod(icount,50).eq.0)) then
          open(UNIT=10)
          call DCHEBYSHEV (u, u1, u2, n)
          do i = 0, n
            eta = etaout*(1-tanh(gamma))*.5*(xi(i)+1)/
     .                   (1-tanh(gamma/2.*(xi(i)+1)))
            write (10,10) eta,u(i),phi1(i)*u1(i),
     .                   (phi1(i)**2*u2(i)+phi2(i)*u1(i))
  10        format (1x,4(e11.5,2x))
          end do
          close(UNIT=10)
        end if
      end do
      write (*,*)
      write (*,*) "Completed solve with interations, residual"
      write (*,*) icount, residual

      return
      end
C******************************************************************************
C******************************************************************************
      SUBROUTINE DCHEBYSHEV(Y, Y1, Y2, N)
C******************************************************************************
C
C     Calculate the Chebyshev transform of a set Y of N+1 real-valued data
C     points and compute the first and second derivatives in Chebyshev space.
C     Then inverse transform and return the derivatives in Y1 and Y2 in real
C     space.  Note that Y is returned unscathed.
C
C******************************************************************************
      REAL Y(0:N), Y1(0:N), Y2(0:N)
      PARAMETER (idim=128)
      REAL WORK(0:IDIM)

      IF (N.GT.IDIM) THEN
        WRITE (*,*) 'ERROR.  N > IDIM in DCHEBYSHEV'
        STOP
      END IF
C
C     SAVE THE INPUT VECTOR
C
      DO I = 0, N
        WORK(I) = Y(I)
      END DO
C
C     COMPUTE THE CHEBYSHEV TRANSFORM
C
      CALL chebyshev(Y,N,1)
C
C     NOW USE THE RECURSIVE RELATION TO TAKE THE FIRST DERIVATIVE
C
      Y1(N) = 0.0
      Y1(N-1) = 2.*FLOAT(N)*Y(N)
      DO K = N-2, 0, -1
        Y1(K) = Y1(K+2) + 2.*(K+1.)*Y(K+1)
      END DO
      Y1(0) = Y1(0)/2.
C
C     NOW REPEAT THE RECURSIVE RELATION TO TAKE THE SECOND DERIVATIVE
C
      Y2(N) = 0.0
      Y2(N-1) = 2.*FLOAT(N)*Y1(N)
      DO K = N-2, 0, -1
        Y2(K) = Y2(K+2) + 2.*(K+1.)*Y1(K+1)
      END DO
      Y2(0) = Y2(0)/2.
C
C     INVERSE TRANSFORM TO GET BACK TO REAL SPACE
C
      CALL chebyshev(Y1,N,-1)
      CALL chebyshev(Y2,N,-1)
      DO I = 0, N
        Y(I) = WORK(I)
      END DO

      RETURN
      END

C***********************************************************************
      subroutine CHEBYSHEV (Y, N, ISIGN)
C***********************************************************************
C
C     Take the Chebyshev transform and normalize
C
C***********************************************************************
      integer  n, isign
      real     y(0:n), t(0:256)

      if (isign .ne. 1) then
        if (isign .ne. -1) then
          write (*,*) 'ERROR:  Invalid Isign in CHEBYSHEV'
          stop
        end if
      end if

      do i = 1, n
        t(i) = y(i)
      end do

      call COSFT3 (y,n,isign)

c     call COSFT (y,n,isign)

      if (isign .eq. 1) then
        do i = 0,n
          y(i) = y(i)/float(n)*2.
        end do
      end if

      return
      end

C******************************************************************************
      SUBROUTINE COSFT3 (Y,N,ISIGN)
C******************************************************************************
C
C     Do a brute force transform.  ISIGN = 1 does a forward transform.
C     ISIGN = -1 a backward transform.  This is a hybrid double/single
C     precision routine.
C
C******************************************************************************
      integer     N, ISIGN
      real        Y(0:N), T(0:N)
      real*16     YY(0:N), TT(0:N), PI

      PI = DACOS(-1.0D0)

      DO I = 0, N
        TT(I) = DBLE(Y(I))
      END DO

      if (isign .eq. 1) then
        DO I = 0, N
          YY(I) = TT(0)/2.D0 + TT(N)*COS(DFLOAT(I)*PI)/2.D0
          DO M = 1, N-1
            YY(I) = YY(I) + TT(M)*COS(DFLOAT(M)*DFLOAT(I)*PI/DFLOAT(N))
          END DO
          YY(I) = YY(I)*2./DFLOAT(N)
        end do
        YY(0) = YY(0)/2.D0
        YY(N) = YY(N)/2.D0
      else
        DO I = 0, N
          YY(I) = 0.0
          DO M = 0, N
            YY(I) = YY(I) + TT(M)*COS(DFLOAT(M)*DFLOAT(I)*PI/DFLOAT(N))
          END DO
        END DO
      end if
C
C     GO BACK TO SINGLE PRECISION
C
      DO I = 0, N
        if (Y(i).le.1e-14) then
          Y(I) = 0.0
        end if
        Y(I) = SNGL(YY(I))
      END DO

      RETURN
      END

C*****************************************************************************
      subroutine TRIDIAG(num,u,a,b,c,f)
C*****************************************************************************
c
c     Solves the tridiagonal system B[a,b,c]u = f
c
c*****************************************************************************
      real    u(num), a(num), b(num), c(num), f(num)

      a(num) = a(num)/b(num)
      do i = num-1,1,-1
        b(i) = b(i)-c(i)*a(i+1)
        a(i) = a(i)/b(i)
      end do

      u(num) = f(num)/b(num)
      do i = num-1,1,-1
        u(i) = (f(i)-c(i)*u(i+1))/b(i)
      end do

      do i = 2, num
        u(i)=u(i)-a(i)*u(i-1)
      end do

      return
      end

C******************************************************************************
      SUBROUTINE ICHEBYSHEV(Y,YI,N)
C******************************************************************************
C
C     Calculate the Chebyshev transform of a set Y of N+1 real-valued data
C     points and integrate in Chebyshev space.  The integral is returned
C     in real space in YI and Y is unscathed.
C
C******************************************************************************
      REAL Y(0:N), YI(0:N+1)
      PARAMETER (idim=128)
      REAL WORK(0:IDIM)

      DO I = 0, N
        WORK(I) = Y(I)
      END DO
C
C     COMPUTE THE CHEBYSHEV TRANSFORM
C
      CALL Chebyshev(Y,N,1)
C
C     NOW INTEGRATE
C     NOTE I AM ASSUMING THAT Y(N) IS SMALL SUCH THAT YI(N+1) = 0.0
C
      YI(N+1) = Y(N)/(2.*(N+1.))
      YI(N) = Y(N-1)/(2.*N)
      YI(1) = 1./2.*(2.*Y(0)-Y(2))
      YI(0) = YI(1)
      DO K = 2, N-1
        YI(K) = 1./(2.*K)*(Y(K-1)-Y(K+1))
        YI(0) = YI(0) + (-1.)**(K-1)*YI(K)
      END DO

      CALL Chebyshev(YI,N,-1)
      DO I = 0, N
        Y(I) = WORK(I)
      END DO

      RETURN
      END
C******************************************************************************
      SUBROUTINE CHEBYINT(N,U,V,NBIG)
C******************************************************************************
C
C
C******************************************************************************
      INTEGER   N, NBIG
      REAL      U(0:N), V(0:N)
      REAL      IU, IV, X, PI

      parameter   (idim=128)
      real        gamma, etaout, phi1(0:idim), phi2(0:idim), beta, alpha
      common      /metrics/  gamma, etaout, phi1, phi2, beta, alpha

      PI = ACOS(-1.0)

      call Chebyshev (u,n,1)
      call Chebyshev (v,n,1)

      open(10)
      do i = 0, NBIG
        X = I*PI/NBIG
        IU = 0.0
        IV = 0.0
        DO M = 0, N
          IU = IU + U(M)*COS(FLOAT(M)*X)
          IV = IV + V(M)*COS(FLOAT(M)*X)
        END DO
        IU = FLOAT(N)*IU*0.5
        IV = FLOAT(N)*IV*0.5
        eta = etaout*(1-tanh(gamma))*.5*(cos(X)+1)/
     .               (1-tanh(gamma/2.*(cos(X)+1)))
        write (10,10) eta, IU, IV
      end do
      close(10)
  10  format (1x,3(e16.8,4x))

      RETURN
      END

C******************************************************************************
      SUBROUTINE SAVE_VELOCITY(N,U,V, xi)
C******************************************************************************
C
C
C******************************************************************************
      INTEGER   N, i, j
      REAL      U(0:N), V(0:N), xi(0:N), eta
      character*15 filename

      parameter   (idim=128)
      real        gamma, etaout, phi1(0:idim), phi2(0:idim), beta, alpha
      common      /metrics/  gamma, etaout, phi1, phi2, beta, alpha

      write (*,5)
   5  format (/,1x,'Write file')
      write (*,10)
  10  format (1x,'Enter filename ==> ',$)
      read (*,'(a)') filename
      open (unit=11,file=filename,status='unknown')
      write (11,*) N
      do i = 0, n
        eta = etaout*(1.-tanh(gamma))*.5*(xi(i)+1.)/
     .               (1.-tanh(gamma/2.*(xi(i)+1.)))
        write (11,*) eta,u(i),v(i)
      end do
      close (11)
      write (*,*)
      write (*,*) 'Wrote file = ', filename

      RETURN
      END
C******************************************************************************
      SUBROUTINE READ_VELOCITY(N,U,V)
C******************************************************************************
C
C
C******************************************************************************
      INTEGER   N, i, j, nmode
      REAL      U(0:N), V(0:N),eta
      character*15 filename

      write (*,5)
   5  format (/,1x,'Read file',/)
      write (*,10)
  10  format (1x,'Enter filename ==> ',$)
      read (*,'(a)') filename
      open (unit=11,file=filename,status='unknown')
      read (11,*) nmode
      do i = 0, n
        read (11,*) eta,u(i),v(i)
      end do
      close (11)

      RETURN
      END

max_stars_repo_name: sscollis/os-stab

max_issues_repo_name: sscollis/os-stab

max_forks_repo_name: sscollis/os-stab


text:
        subroutine e3mtrx (rho,    pres,   T,
     &                     ei,     h,	   alfap,
     &                     betaT,  cp,     rk,
     &                     u1,     u2,   u3,
     &                     um1,    um2,  um3,
     &                     A0,     A1,
     &                     A2,     A3,
     &                     e2p,    e3p,    e4p,
     &                     drdp,   drdT,   A0DC,
     &                     A0inv,  dVdY)
c
c----------------------------------------------------------------------
c
c This routine sets up the necessary matrices at the integration point.
c
c input:
c  rho   (npro)         : density
c  pres  (npro)         : pressure
c  T     (npro)         : temperature
c  ei    (npro)         : internal energy
c  h     (npro)         : enthalpy
c  alfap (npro)         : expansivity
c  betaT (npro)         : isothermal compressibility
c  cp    (npro)         : specific heat at constant pressure
c  c     (npro)         : speed of sound
c  rk    (npro)         : kinetic energy
c  u1    (npro)         : x1-velocity component
c  u2    (npro)         : x2-velocity component
c  u3    (npro)         : x3-velocity component
c  um1    (npro)         : ALE x1-velocity component
c  um2    (npro)         : ALE x2-velocity component
c  um3    (npro)         : ALE x3-velocity component
c
c output:
c  A0    (npro,nflow,nflow)  : A0 matrix
c  A1   (npro,nflow,nflow)  : A_1 matrix
c  A2   (npro,nflow,nflow)  : A_2 matrix
c  A3   (npro,nflow,nflow)  : A_3 matrix
c
c Note: the definition of the matrices can be found in
c       thesis by Hauke.
c
c Zdenek Johan, Summer 1990.  (Modified from e2mtrx.f)
c Zdenek Johan, Winter 1991.  (Fortran 90)
c Kenneth Jansen, Winter 1997 Primitive Variables
c----------------------------------------------------------------------
c
        include "common.h"
c
c
c  passed arrays
c
        dimension rho(npro),                 pres(npro),
     &            T(npro),                   ei(npro),
     &            h(npro),                   alfap(npro),
     &            betaT(npro),
     &            cp(npro),
     &            rk(npro),
     &            u1(npro),                 u2(npro),
     &            u3(npro),                 fact1(npro),
     &            um1(npro), um2(npro), um3(npro),
     &            A0(npro,nflow,nflow),     dVdY(npro,15),
     &            A1(npro,nflow,nflow),     A2(npro,nflow,nflow),
     &            A3(npro,nflow,nflow),     A0DC(npro,4),
     &            A0inv(npro,15),           d(npro),
     &            fact2(npro),               s1(npro),
     &            e1bar(npro),              e2bar(npro),
     &            e3bar(npro),              e4bar(npro),
     &            e5bar(npro),              c1bar(npro),
     &            c2bar(npro),              cv(npro),
     &            c3bar(npro),              u12(npro),
     &            u31(npro),                u23(npro)
c
c  local work arrays that are passed shared space
c
        dimension e2p(npro),
     &            e3p(npro),                 e4p(npro),
     &            drdp(npro),                drdT(npro)

	ttim(21) = ttim(21) - secs(0.0)
c
c.... initialize
c
        A0 = zero
        A1 = zero
        A2 = zero
        A3 = zero
c
c.... set up the constants
c
c
        drdp = rho * betaT
        drdT = -rho * alfap
        A0(:,5,1) = drdp * (h + rk)  - alfap * T    ! e1p
c        A0(:,5,1) = drdp * (ei + rk) + betaT * pres - alfap * T    ! e1p
          e2p  = A0(:,5,1) + one
          e3p  = rho * ( h + rk)
          e4p  = drdT * (h + rk) + rho * cp
c
c
c.... Calculate A0
c
        A0(:,1,1) = drdp
c       A0(:,1,2) = zero
c       A0(:,1,3) = zero
c       A0(:,1,4) = zero
        A0(:,1,5) = drdT
c
        A0(:,2,1) = drdp * u1
        A0(:,2,2) = rho
c       A0(:,2,3) = zero
c       A0(:,2,4) = zero
        A0(:,2,5) = drdT * u1
c
        A0(:,3,1) = drdp * u2
c       A0(:,3,2) = zero
        A0(:,3,3) = rho
c       A0(:,3,4) = zero
        A0(:,3,5) = drdT * u2
c
        A0(:,4,1) = drdp * u3
c       A0(:,4,2) = zero
c       A0(:,4,3) = zero
        A0(:,4,4) = rho
        A0(:,4,5) = drdT * u3
c
covered above       A0(:,5,1) = drdp * u1
        A0(:,5,2) = rho * u1
        A0(:,5,3) = rho * u2
        A0(:,5,4) = rho * u3
        A0(:,5,5) = e4p
c
   !      flops = flops + 67*npro
c
c.... ALE version of A_1, A_2 and A_3
c.... Calculate A-tilde-1, A-tilde-2 and A-tilde-3
c
        A1(:,1,1) = drdp * (u1 - um1)
        A1(:,1,2) = rho
c       A1(:,1,3) = zero
c       A1(:,1,4) = zero
        A1(:,1,5) = drdT * (u1 - um1)
c
        A1(:,2,1) = drdp * (u1 - um1) * u1 +1
        A1(:,2,2) = rho  * (u1 - um1) + rho * u1
c       A1(:,2,3) = zero
c       A1(:,2,4) = zero
        A1(:,2,5) = drdT * (u1 - um1) * u1
c
        A1(:,3,1) = drdp * (u1 - um1) * u2
        A1(:,3,2) = rho  * u2
        A1(:,3,3) = rho  * (u1 - um1)
c       A1(:,3,4) = zero
        A1(:,3,5) = drdT * (u1 - um1) * u2
c
        A1(:,4,1) = drdp * (u1 - um1) * u3
        A1(:,4,2) = rho  * u3
c       A1(:,4,3) = zero
        A1(:,4,4) = rho  * (u1 - um1)
        A1(:,4,5) = drdT * (u1 - um1) * u3
c
c        A1(:,5,1) = (u1 - um1) * e2p
      A1(:,5,1) = u1*e2p - um1*A0(:,5,1)
        A1(:,5,2) = e3p + rho * (u1 - um1) * u1
        A1(:,5,3) = rho * (u1 - um1) * u2
        A1(:,5,4) = rho * (u1 - um1) * u3
        A1(:,5,5) = (u1 - um1) * e4p
c
   !      flops = flops + 35*npro
c
        A2(:,1,1) = drdp * (u2 - um2)
c       A2(:,1,2) = zero
        A2(:,1,3) = rho
c       A2(:,1,4) = zero
        A2(:,1,5) = drdT * (u2 - um2)
c
        A2(:,2,1) = drdp * u1 * (u2 - um2)
        A2(:,2,2) = rho  * (u2 - um2)
        A2(:,2,3) = rho  * u1
c       A2(:,2,4) = zero
        A2(:,2,5) = drdT * u1 * (u2 - um2)
c
        A2(:,3,1) = drdp * (u2 - um2) * u2 + 1
c       A2(:,3,2) = zero
        A2(:,3,3) = rho  * (u2 - um2) + rho * u2
c       A2(:,3,4) = zero
        A2(:,3,5) = drdT * (u2 - um2) * u2
c
        A2(:,4,1) = drdp * (u2 - um2) * u3
c       A2(:,4,2) = zero
        A2(:,4,3) = rho  * u3
        A2(:,4,4) = rho  * (u2 - um2)
        A2(:,4,5) = drdT * (u2 - um2) * u3
c
c        A2(:,5,1) = (u2 - um2) * e2p
        A2(:,5,1) = u2 * e2p - um2 * A0(:,5,1)
        A2(:,5,2) = rho * u1 * (u2 - um2)
        A2(:,5,3) = e3p + rho * (u2 - um2) * u2
        A2(:,5,4) = rho * (u2 - um2) * u3
        A2(:,5,5) = (u2 - um2) * e4p
c
   !      flops = flops + 35*npro
c
        A3(:,1,1) = drdp * (u3 - um3)
c       A3(:,1,2) = zero
c       A3(:,1,3) = zero
        A3(:,1,4) = rho
        A3(:,1,5) = drdT * (u3 - um3)
c
        A3(:,2,1) = drdp * u1 * (u3 - um3)
        A3(:,2,2) = rho  * (u3 - um3)
c       A3(:,2,3) = zero
        A3(:,2,4) = rho  * u1
        A3(:,2,5) = drdT * u1 * (u3 - um3)
c
        A3(:,3,1) = drdp * (u3 - um3) * u2
c       A3(:,3,2) = zero
        A3(:,3,3) = rho  * (u3 - um3)
        A3(:,3,4) = rho  * u2
        A3(:,3,5) = drdT * (u3 - um3) * u2
c
        A3(:,4,1) = drdp * (u3 - um3) * u3 + 1
c       A3(:,4,2) = zero
c       A3(:,4,3) = zero
        A3(:,4,4) = rho  * u3 + rho * (u3 - um3)
        A3(:,4,5) = drdT * (u3 - um3) * u3
c
c        A3(:,5,1) = (u3 - um3) * e2p
        A3(:,5,1) = u3 * e2p - um3 * A0(:,5,1)
        A3(:,5,2) = rho * u1 * (u3 - um3)
        A3(:,5,3) = rho * u2 * (u3 - um3)
        A3(:,5,4) = e3p + rho * (u3 - um3) * u3
        A3(:,5,5) = (u3 - um3) * e4p
c
   !      flops = flops + 35*npro
	ttim(21) = ttim(21) + secs(0.0)

c
c.... return
c
      if (idc .ne. 0) then
c.... for Discountinuity Capturing Term
c
c.... calculation of A0^DC matrix
c
c.... Ref P-163 of the handout
c
       s1 = one/(rho**2 * betaT * T)
       cv = cp - (alfap**2 * T/rho/betaT)
       A0DC(:,1) = (rho * betaT)**2*s1
       A0DC(:,2) = -rho*alfap*rho*betaT*s1
       A0DC(:,3) = rho/T
       A0DC(:,4) = (-rho*alfap)**2 * s1 + (rho*cv/T**2)
c
c.... calculation of A0^tilda^inverse matrix
c.... ref P-169 of the hand out
c
       fact1 = one/(rho*cv*T**2)
       d = alfap*T/rho/betaT
       e1bar = h - rk
       e2bar = e1bar - d
       e3bar = e2bar - cv * T
       e4bar = e2bar - 2* cv * T
       e5bar = e1bar**2 - 2*e1bar*d + 2*rk*cv*T + cp*T/rho/betaT
       c1bar = u1**2 + cv * T
       c2bar = u2**2 + cv * T
       c3bar = u3**2 + cv * T
       u12 = u1 * u2
       u31 = u3 * u1
       u23 = u2 * u3
       A0inv( :,1) = e5bar*fact1
       A0inv( :,2) = c1bar*fact1
       A0inv( :,3) = c2bar*fact1
       A0inv( :,4) = c3bar*fact1
       A0inv( :,5) = 1*fact1
       A0inv( :,6) = u1*e3bar*fact1
       A0inv( :,7) = u2*e3bar*fact1
       A0inv( :,8) = u3*e3bar*fact1
       A0inv( :,9) = -e2bar*fact1
       A0inv(:,10) = u12*fact1
       A0inv(:,11) = u31*fact1
       A0inv(:,12) = -u1*fact1
       A0inv(:,13) = u23*fact1
       A0inv(:,14) = -u2*fact1
       A0inv(:,15) = -u3*fact1
c
c.....calculation of dV/dY (derivative of entropy variables w.r.to primitive
c
      fact1 = 1/T
      fact2 = fact1/T
      dVdY( :,1) = fact1/rho
      dVdY( :,2) = -fact1*u1
      dVdY( :,3) =  fact1
      dVdY( :,4) =  -fact1*u2
      dVdY( :,5) =  zero
      dVdY( :,6) =  fact1
      dVdY( :,7) =  -fact1*u3
      dVdY( :,8) =  zero
      dVdY( :,9) =  zero
      dVdY(:,10) =  fact1
      dVdY(:,11) =  -(h-rk)*fact2
      dVdY(:,12) =  -fact2*u1
      dVdY(:,13) =  -fact2*u2
      dVdY(:,14) =  -fact2*u3
      dVdY(:,15) =   fact2

      endif  !end of idc.ne.0

        return
        end
c
c
        subroutine e3mtrxSclr (rho,
     &                         u1,     u2,   u3,
     &                         A0t,    A1t,
     &                         A2t,    A3t  )
c
c----------------------------------------------------------------------
c
c This routine sets up the necessary matrices at the integration point.
c
c input:
c  rho   (npro)         : density
c  u1    (npro)         : x1-velocity component
c  u2    (npro)         : x2-velocity component
c  u3    (npro)         : x3-velocity component
c
c output:
c  A0t    (npro) : A_0 "matrix"
c  A1t   (npro)  : A_1 "matrix"
c  A2t   (npro)  : A_2 "matrix"
c  A3t   (npro)  : A_3 "matrix"
c
c Note: the definition of the matrices can be found in
c       thesis by Hauke.
c
c Zdenek Johan, Summer 1990.  (Modified from e2mtrx.f)
c Zdenek Johan, Winter 1991.  (Fortran 90)
c Kenneth Jansen, Winter 1997 Primitive Variables
c----------------------------------------------------------------------
c
        include "common.h"
c
c
c  passed arrays
c
        dimension rho(npro),
     &            u1(npro),        u2(npro),
     &            u3(npro),
     &            A0t(npro),
     &            A1t(npro),       A2t(npro),
     &            A3t(npro)
c
        if (iconvsclr.eq.2) then  !convective form
           A0t(:) = one
           A1t(:) = u1(:)
           A2t(:) = u2(:)
           A3t(:) = u3(:)
        else                    !conservative form
           A0t(:) = rho(:)
           A1t(:) = rho(:)*u1(:)
           A2t(:) = rho(:)*u2(:)
           A3t(:) = rho(:)*u3(:)
        endif

c
c.... return
c
        return
        end


text:
!	print eap data for latent vector.
!	comment is used for identification.
!	factorname gives names of elements of the latent vector.
!	id is individual id.
!	unit uniteap is the unit for output.
!	eapid indicates if ids are put out.
!	covtheta are conditional covariance matrices.
!	meantheta are eaps.

subroutine eapoutput(comment,factorname,id,uniteap,eapid,covtheta,meantheta)
implicit none
character(len=*),intent(in)::comment
character(len=32),intent(in)::factorname(:)
character(len=32),intent(in),optional::id(:)
integer,intent(in)::uniteap
logical,intent(in)::eapid
real(kind=8),intent(in)::covtheta(:,:,:),meantheta(:,:)

!	Buffers for writing.
character(len=4)::buffer
character(len=32)::bufferid
character(len=25),allocatable::buffm(:),buffc(:,:)
!	format
character(len=7)::writefmt
!	col is column number.
!	io is error flag.
!	obs is observation number.
!	row is row number.
integer::col,io,obs,row
col=size(factorname)

allocate(buffm(col),buffc(col,col),stat=io)
if(io/=0) stop "Allocation failure for arrays for output of conditional means and covariances."
write(uniteap,'(2a)') "Conditional expectations and covariances ",comment
write(buffer,'(i4)') 1+2*col*(col+1)
writefmt='('//trim(adjustl(buffer))//'a)'
write(unit=uniteap,fmt=writefmt,iostat=io) "ID",(",","Mean_"//trim(adjustl(factorname(row))),row=1,size(factorname)),&
	((",","Cov_"//trim(adjustl(factorname(row)))//"_"//trim(adjustl(factorname(col))),&
	col=1,size(factorname)),row=1,size(factorname))
if(io/=0) stop "Writing of conditional means and covariances failed."

do obs=1,size(meantheta,2)
	do row=1,size(factorname)
		write(buffm(row),'(g25.16e3)')meantheta(row,obs)
		do col=1,size(factorname)
			write(buffc(row,col),'(g25.16e3)')covtheta(row,col,obs)
		end do
	end do
	if(eapid)then
		bufferid=id(obs)
	else
		write(bufferid,'(i12)') obs
	end if
	write(unit=uniteap,fmt=writefmt,iostat=io) trim(adjustl(bufferid)),(",",trim(adjustl(buffm(row))),row=1,size(factorname)),&
			((",",trim(adjustl(buffc(row,col))),col=1,size(factorname)),row=1,size(factorname))

	if(io/=0) stop "Writing of conditional means and covariances failed."
end do
return

end subroutine eapoutput

max_stars_repo_name: EducationalTestingService/MIRT

max_issues_repo_name: EducationalTestingService/MIRT

max_forks_repo_name: EducationalTestingService/MIRT


text:
      subroutine renorm (dexv, vcoul, srho)

      implicit double precision (a-h,o-z)
      save

      common /print/ iprint
      common /atomco/ den(30), dq1(30), dfl(30), ws, nqn(30), nql(30),
     1                nk(30), nmax(30), nel(30), norb, norbco
      common /dira/ dv(251), dr(251), dp(251), dq(251), dpas, tets,
     1              z, nstop, nes, np, nuc
      common /deux/ dvn(251), dvf(251), d(251), dc(251), dgc(251,30),
     1 dpc(251,30)

c     vcoul is the coulomb potential (no factor of r**2) (output)
      dimension vcoul(251)
c     srho is charge density in form 4*pi*density*r**2 output)
      dimension srho(251)
c jm  9/23/87 added srho renormalized charge density to be used
c     in cphase
       external exchan
      do 10 i=1,np
         dv(i)=0.0
         d(i)=0.0
   10 continue
      ddjri=log(ws/dr(1))/dpas
      jri=1.0+ddjri
      jr1=jri
      ddjr1=ddjri-jr1+1.0

      if (jri-2*(jri/2).ne.0) go to 20
         jri=jri+1
   20 continue

      ddjri=ddjri-jri+1.0
c  ddjri = (log(ws)-dri)/dpas
c  dri  =  log(dr(jri))

      da=0.0
      do 30 j=1,norb
      do 30 i=1,np
   30    d(i)=d(i)+nel(j)*(dgc(i,j)**2+dpc(i,j)**2)

      do 50 i=jri,np
         dl=dr(i)
         if (i.eq.jri.or.i.eq.np) go to 40
            dl=dl+dl
            if ((i-2*(i/2)).eq.0) dl=dl+dl
   40    dd=d(i)*dl
         da=da+dd
   50 continue

      da=dpas*da/3.0
      dfo=dr(jri-1)*d(jri-1)
      df1=dr(jri)*d(jri)
      df2=dr(jri+1)*d(jri+1)
      dcor=-dpas*(df1*ddjri+(df2+dfo-2.0*df1)*ddjri**3/6.0+(df2-dfo)
     1 *ddjri**2*.25)
      da=da+dcor
      if (iprint .ge. 5)  write(16,60) da
   60 format (1h ,' no. of electrons outside the ws-radius',e16.8)
      db=0.0

      do 80 i=jri,np
         dl=1.0
         if (i.eq.jri.or.i.eq.np) go to 70
            dl=dl+dl
            if ((i-2*(i/2)).eq.0) dl=dl+dl
   70    dd=d(i)*dl
         db=db+dd
   80 continue

      db=dpas*db/3.0
      df0=d(jri-1)
      df1=d(jri)
      df2=d(jri+1)
      dcor=-dpas*(df1*ddjri+(df2+df0-2.0*df1)*ddjri**3/6.0+(df2-df0)
     1 *ddjri**2*.25)
      db=db+dcor
      if (iprint .ge. 5)  write(16,90) db
   90 format (1h ,' db= ',e16.8)

      call potslw (dvn,d,dp,dr,dpas,np)

      du=da*3.0/(ws**3)

      do 120 i=1,np
         if (i.gt.jr1+1) then
            srho(i)=0.0
            go to 100
         endif
            d(i)=d(i)+du*dr(i)**2
            srho(i)=d(i)
  100    continue
         dumm=-exchan(d(i),dr(i),dexv)/dr(i)
         dvf(i)=dumm
         if (i.gt.jr1) go to 110
            dvn(i)=dvn(i)-z/dr(i)+da*(1.50/ws-.50*dr(i)**2/ws**3)-db
            go to 120
  110    continue
            dvn(i)=0.0
  120 dv(i)=dvn(i)+dumm

c ad1 write the mt index and radius
      if (iprint .ge. 5)  write(16,55)jr1,dr(jr1)
  55  format(' jr1 = ',i10,10x,'wigner-seitz radius = ',e16.8)

c ad1 output 2.*dvn*r**2 for use in phase (dvn = normalised coulomb)
c     write(17,200)((2.0*dvn(i)*dr(i)*dr(i)),i=1,np)
c 200 format(1p5e16.8)
c      passvc formerly used to pass data directly to PHASE
c      do 151  i = 1, np
c         passvc (i) = 2.0 * dvn(i) * dr(i) * dr(i)
c  151 continue
c
c     passvc above is vcoul*r**2
      do 151  i = 1, np
         vcoul(i) = 2 * dvn(i)
  151 continue

c jm  output renormalized charge density for use in cphase
c                                          (d=4pi*rho*r^2)
c     write(18,200) srho

cjm write out rs as function of r
c     do 8934 i=1,jr1
c     xxrs=(3*dr(i)*dr(i)/srho(i))**.33333333
c8934 write(19,140) dr(i), xxrs
      return
      end

max_stars_repo_name: keechul/ifeffit

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: bruceravel/ifeffit

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: bruceravel/ifeffit

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
module Wave_Packet_Simulation

    use iso_fortran_env, only : REAL64
    use Wavefunc_DataStorage
    use omp_lib
    implicit none

    contains

    ! Simple wave packet centered at (0,0) with momentum px, py and width sigma
    pure function start_wave_packet(x, y, px, py, sigma, hbar, pi) result (psi_at_point)
        REAL(kind=REAL64), intent(in) :: x, y, px, py, sigma, hbar, pi
        COMPLEX(kind=REAL64) :: psi_at_point

        REAL(kind=REAL64) :: coeff_val
        COMPLEX(kind=REAL64) :: exp_val

        !coeff_val = sqrt(sigma *  sqrt(pi))
        !exp_val = (px * x + py * y) * (0, 1.0)/hbar - (x * x + y * y)/(2 * sigma * sigma)

        !psi_at_point = coeff_val * exp(exp_val)

        if (x .eq. 0 .and. y .eq. 0) then
            psi_at_point = 1
        else
            psi_at_point = 0
        end if

    end function start_wave_packet

    ! The first step iteration, since we don't have access to t-s at that point
    pure function first_iteration_at_point(xind, yind, s, d, hbar, m, curV, curPsi) result(psi_next)
        Integer, intent(in) :: xind, yind
        Real(kind=REAL64), intent(in) :: s, d, hbar, m, curV
        Complex(kind=REAL64), intent(in), dimension(:, :) :: curPsi

        Complex(kind=REAL64) :: psi_next

        Complex(kind=REAL64) :: offPsi, offPsi_coeff, psiV_coeff, psi_coeff, curPsiPoint
        Integer, dimension(2) :: realitySize

        psi_coeff = 1 - (0, -2) *hbar * s/(m * d * d)
        psiV_coeff = (0, -1) * s /hbar
        offPsi_coeff = (0, 0.5) * hbar * s/(m * d * d)

        curPsiPoint = curPsi(xind, yind)

        ! We default the value of Psi to 0 at the borders, since we are assuming
        ! a particle in a finite box.
        realitySize = shape(curPsi)
        if (xind .gt. 1 .and. xind .lt. realitySize(1)) then
            offPsi = curPsi(xind + 1, yind) + curPsi(xind - 1, yind)
        else if(xind .lt. realitySize(1)) then
            offPsi = curPsi(xind + 1, yind) + 0
        else
            offPsi = curPsi(xind - 1, yind) + 0
        end if

        ! Same for y
        if (yind .gt. 1 .and. yind .lt. realitySize(2)) then
            offPsi = offPsi + curPsi(xind, yind + 1) + curPsi(xind, yind - 1)
        else if(yind .lt. realitySize(2)) then
            offPsi = offPsi + curPsi(xind, yind + 1) + 0
        else
            offPsi = offPsi + curPsi(xind, yind - 1) + 0
        end if

        psi_next = curPsiPoint * psi_coeff + curPsiPoint * curV * psiV_coeff + offPsi * offPsi_coeff

    end function first_iteration_at_point

    ! Get the coefficient for the \Psi(x, y, t) term
    pure function calc_psi_coeff(m, hbar, s, d) result(retval)
        REAL(kind=REAL64), intent(in) :: m, hbar, s, d
        COMPLEX(kind=REAL64) :: retval

        retval = (-4.0, 1.0) * s * hbar / (m * d * d)

    end function calc_psi_coeff

    ! Get the coefficient for the \V(x, y, t) * \Psi(x, y, t) term
    pure function calc_psiV_coeff(m, hbar, s, d) result(retval)
        REAL(kind=REAL64), intent(in) :: m, hbar, s, d
        COMPLEX(kind=REAL64) :: retval

        retval = (-2.0, 1.0) * s/hbar

    end function calc_psiV_coeff

    ! Get the coefficient for the gradient components summation term
    pure function calc_psiOff_coeff(m, hbar, s, d) result(retval)
        REAL(kind=REAL64), intent(in) :: m, hbar, s, d
        COMPLEX(kind=REAL64) :: retval

        retval = hbar * s * (0.0, 1.0) / (m * d * d)
    end function calc_psiOff_coeff

    ! General function to iterate the wave at a point
    pure function iterate_psi_at_point(xind, yind, psi_coeff, psiV_coeff, psiOff_coeff, curV, psi_prev, psi_cur) result(psi_next)
        INTEGER, intent(in) :: xind, yind
        REAL(kind=REAL64), intent(in) :: curV
        COMPLEX(kind=REAL64), intent(in) :: psi_coeff, psiV_coeff, psiOff_coeff
        Complex(kind=REAL64), dimension(:, :), intent(in) :: psi_prev, psi_cur
        COMPLEX(kind=REAL64) :: psi_next

        COMPLEX(kind=REAL64) :: offPsi
        Integer, dimension(2) :: realitySize

        ! We default the value of Psi to 0 at the borders, since we are assuming
        ! a particle in a finite box.
        realitySize = shape(psi_cur)
        if (xind .gt. 1 .and. xind .lt. realitySize(1)) then
            offPsi = psi_cur(xind + 1, yind) + psi_cur(xind - 1, yind)
        else if(xind .lt. realitySize(1)) then
            offPsi = psi_cur(xind + 1, yind) + 0
        else
            offPsi = psi_cur(xind - 1, yind) + 0
        end if

        ! Same for y
        if (yind .gt. 1 .and. yind .lt. realitySize(2)) then
            offPsi = offPsi + psi_cur(xind, yind + 1) + psi_cur(xind, yind - 1)
        else if(yind .lt. realitySize(2)) then
            offPsi = offPsi + psi_cur(xind, yind + 1) + 0
        else
            offPsi = offPsi + psi_cur(xind, yind - 1) + 0
        end if

        psi_next = psi_prev(xind, yind) + (psi_coeff + curV * psiV_coeff) * psi_cur(xind, yind) + psiOff_coeff * offPsi

    end function iterate_psi_at_point

    ! Normalizes the wave function data, so sum(abs(matrix)**2) = 1.
    subroutine normalize_matrix(matrix)
        COMPLEX(kind=REAL64), intent(inout), dimension(:, :) :: matrix

        INTEGER, dimension(2) :: dims
        REAL(kind=REAL64) :: sum_squares
        INTEGER :: xind, yind

        dims = SHAPE(matrix)
        sum_squares = 0.0

        ! Get the normalization coefficient at that time
        ! $OMP parallel do default(shared) private(xind, yind) reduction (+:sum_squares)
        do xind = 1, dims(1)
            do concurrent (yind = 1:dims(2))
                if (matrix(xind, yind) .ne. 0) then
                    sum_squares = sum_squares + 1D100 * custom_abs(matrix(xind, yind))**2
                end if
            end do
        end do
        ! $OMP end parallel do

        sum_squares = sum_squares /1D100

        ! Short Circuit on a bad universe
        if (sum_squares .eq. 0) then
            Print *, "Bad universe. "
            return
        else if(abs(sum_squares - 1) .lt. 0.0001) then
            Print *, "SPICY! Sumsq = ", sum_squares
        end if

        ! Get the normalization coefficient at that time

        ! $OMP parallel do default(shared) private(xind, yind) shared(sum_squares, matrix)
        do xind = 1, dims(1)
            do concurrent (yind = 1:dims(2))
                if (matrix(xind, yind) .ne. 0) then
                    matrix(xind, yind) = matrix(xind, yind) / sqrt(sum_squares)
                end if
            end do
        end do
        ! $OMP end parallel do

    end subroutine

end module Wave_Packet_Simulation

program main

    use Wave_Packet_Simulation

    implicit none

    ! Constants of reality
    REAL(kind=REAL64), parameter :: PI  = 4.D0 * DATAN(1.D0) ! Gurantees max accuracy
    REAL(kind=REAL64), parameter :: hbar = 1.054571800D-34!J * s

    ! Constants of this simulation

    ! Data storage variables
    type(DataStore) ::  storage
    Integer, parameter :: timesize = 200, points_per_time = 100000, storage_sz = points_per_time * timesize
    Integer, parameter :: max_points = 400000
    Integer, parameter :: xsize = 1000, ysize = 1000 ! Universe size
    COMPLEX(kind=REAL64), dimension(:, :), allocatable, target :: universe_a, universe_b, universe_c ! The closest time slices
    COMPLEX(kind=REAL64), dimension(:, :), pointer :: prev_data, curr_data, next_data, tmp_ptr ! Pointers to each slice to more efficiently rotate

    ! Universe bounds
    REAL(kind=REAL64), parameter :: xmin = -1.0, xmax = 1.0, ymin = -1.0, ymax = 1.0
    REAL(kind=REAL64), parameter :: tmax = 1.0
    REAL(kind=REAL64), parameter :: min_percent = 1D-6

    REAL(kind=REAL64), parameter :: s = (tmax)/timesize ! time step, seconds
    REAL(kind=REAL64), parameter :: d = MAX((xmax - xmin)/xsize, (ymax -ymin)/ysize) ! distance step, meters; equal to equal to (xmax - xmin)/xsize

    REAL(kind=REAL64), parameter :: px = 1.0, py = 0 ! initial momentum, m/s
    REAL(kind=REAL64), parameter :: sigma = 4D-4 ! the width in meters?

    REAL(kind=REAL64), parameter :: m = 9.10938356D-31 ! Mass of particle (electron here) in kg

    Integer, parameter :: fd = 1 ! Output CSV file descriptor

    ! Local variables
    Integer :: xind, yind, tind ! Iteration indices
    REAL(kind=REAL64) :: x, y, t ! Temp variables for converting from the index
    REAL(kind=REAL64) :: curV ! Temp variable for the potential at a point
    COMPLEX(kind=REAL64) :: psi_coeff, psiV_coeff, psi_off_coeff ! Coefficients for the iterator
    type(TimeNode), pointer :: data_idx

    storage%size = 0
    Nullify(storage%head)
    Nullify(storage%tail)

    allocate(universe_a(xsize, ysize))
    allocate(universe_b(xsize, ysize))
    allocate(universe_c(xsize, ysize))

    Print *, "Hello!"
    Print *, "Going from t=0 to t=", timesize * s, " with steps of ", s, &
             " over a region of [",xmin, ", ", xmin + d*xsize, "] x [", ymin, ", ", ymin + d*ysize, &
             "] with steps of ", d

    tmp_ptr => universe_a
    prev_data => universe_a ! The universe at time t+s
    curr_data => universe_b ! The universe at time t
    next_data => universe_c ! The universe at time t+s

    !$omp parallel do default(shared) private(x, y, xind, yind)
    do xind = 1,xsize
        x = (xind - 1) * d + xmin
        do concurrent (yind = 1:ysize)
            y = (yind - 1) * d + ymin
            prev_data(xind, yind) = start_wave_packet(x, y, px, py, sigma, hbar, pi)
        end do
    end do
    !$omp end parallel do

    call normalize_matrix(prev_data)
    call store_data(storage, prev_data, 0D0, d, xmin, ymin, max_points, min_percent)

    Print *, "Finished t = ", 0, " step with ", storage%size, " total points."

    !$omp parallel do default(shared) private(x, y, xind, yind, curV)
    do xind = 1,xsize
        x = (xind - 1) * d + xmin
        do concurrent (yind = 1:ysize)
            y = (yind - 1) * d + ymin
            curV = V(x, y, s)
            curr_data(xind, yind) = first_iteration_at_point(xind, yind, s, d, hbar, m, curV, prev_data)
        end do
    end do
    !$omp end parallel do

    call normalize_matrix(curr_data)

    call store_data(storage, curr_data, s, d, xmin, ymin, max_points, min_percent)

    Print *, "Finished t = ", 1, " step with ", storage%size, " total points."

    psi_coeff = calc_psi_coeff(m, hbar, s, d)
    psiV_coeff = calc_psiV_coeff(m, hbar, s, d)
    psi_off_coeff = calc_psiOff_coeff(m, hbar, s, d)

    tind = 3
    t = s
    do while(storage%size .le. storage_sz .and. t .lt. tmax)
        t = (tind - 1) * s

        !$omp parallel do default(shared) private(x, y, xind, yind, curV)
        do xind = 1,xsize
            x = (xind - 1) * d + xmin
            do concurrent (yind = 1:ysize)
                y = (yind - 1) * d + ymin
                curV = V(x, y, t)
                next_data(xind, yind) = iterate_psi_at_point( &
                    xind, yind, psi_coeff, psiV_coeff, psi_off_coeff, curV,  &
                    prev_data, curr_data &
                )
            end do
        end do
        !$omp end parallel do

        call normalize_matrix(next_data)
        call store_data(storage, next_data, t, d, xmin, ymin, max_points, min_percent)
        tmp_ptr => prev_data
        prev_data => curr_data
        curr_data => next_data
        next_data => tmp_ptr
        if(MOD(tind-1, 100) .eq. 0) then
            Print *, "Finished t = ", tind-1, " step with ", storage%size, " total data points."
        end if
        tind = tind + 1
    end do

    ! Write the data to a disk
    ! To save space we only write the non-zero data, not the entire
    ! universe at every time point

    open(fd, file="Simulationdata.csv")

    data_idx => storage%head
    do while(ASSOCIATED(data_idx))

        do xind = 1, size(data_idx%data)
            Write (fd, *) data_idx%t, ", ", &
                        data_idx%data(xind)%x, ", ", &
                        data_idx%data(xind)%y, ", ", &
                        real(data_idx%data(xind)%value), ",", &
                        aimag(data_idx%data(xind)%value)
        end do
        data_idx => data_idx%next
    end do

    deallocate(universe_a)
    deallocate(universe_b)
    deallocate(universe_c)

    contains

    ! The potential function
    pure function V(x, y, t) result(retval)
        REAL(kind=REAL64), intent(in) :: x, y, t
        REAL(kind=REAL64) :: retval

        ! BIG potential at the edges of the box
        if (t .ge. 0 .and. (x .gt. 1 .or. y .gt. 1 .or. x .lt. -1 .or. y .lt. -1)) then
            retval = 9D99
        else

            ! Right now we have no potential inside
            retval = 0
        end if

    end function V

end program main

max_stars_repo_name: ischeinkman/Quantum-Experiments

max_issues_repo_name: ischeinkman/Quantum-Experiments

max_forks_repo_name: ischeinkman/Quantum-Experiments