text:
       subroutine sum_paths(idata, iupath, nxpath, nqmax,
     $      ckreal, ckimag)
c
c//////////////////////////////////////////////////////////////////////
c Copyright (c) 1997--2000 Matthew Newville, The University of Chicago
c Copyright (c) 1992--1996 Matthew Newville, University of Washington
c
c Permission to use and redistribute the source code or binary forms of
c this software and its documentation, with or without modification is
c hereby granted provided that the above notice of copyright, these
c terms of use, and the disclaimer of warranty below appear in the
c source code and documentation, and that none of the names of The
c University of Chicago, The University of Washington, or the authors
c appear in advertising or endorsement of works derived from this
c software without specific prior written permission from all parties.
c
c THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
c EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
c MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
c IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
c CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
c TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
c SOFTWARE OR THE USE OR OTHER DEALINGS IN THIS SOFTWARE.
c//////////////////////////////////////////////////////////////////////
c
c  sum-over-paths for a particular data set, returning
c      real and imaginary parts of chi(k)
c     used by iff_ff2chi and fitfun
c
       implicit none
       include 'consts.h'
       integer       idata, i, ip, nqmax, xafs_path, ret
       integer       nxpath, iupath(mpaths)
       double precision  xdsave, xtmp, getsca, reff
       double precision  ckreal(maxpts), ckimag(maxpts)
       double precision  tmpi(maxpts), tmpr(maxpts)
       external    getsca, xafs_path

       do 100 i = 1, maxpts
          ckreal(i)  = zero
          ckimag(i)  = zero
 100   continue
       xdsave = getsca('data_set',0)
       xtmp   = one*max(1,min(mdata,idata))
       call setsca('data_set',  xtmp)
cc       print*, 'sum paths data_set', xtmp
       do 1000 ip = 1, nxpath
          ret = xafs_path(iupath(ip), tmpr, tmpi, reff)
cc          print*, ' path ', ip, iupath(ip), ret, reff
          if (ret .eq. 1) then
c  add this to the other paths
             do 850 i = 1, nqmax
                ckreal(i) = ckreal(i) + tmpr(i)
                ckimag(i) = ckimag(i) + tmpi(i)
 850         continue
c          else
c             print*, 'sum over paths: path not used: ',
c     $            ip, iupath(ip)
          endif
 1000  continue
c now return previous the previous 'user variables'
       call setsca('data_set',  xdsave)
       return
       end

       integer function xafs_path(ipath,chi_r, chi_i,reff)
c
c  calculate xafs for a single path
c  ipath =  'user path index'
c  idata =  'data set'
c  reff <= 0. on output means path isn't defined.
c
       implicit none
       include 'consts.h'
       include 'keywrd.h'
       include 'arrays.h'
       include 'fefdat.h'
       include 'feffit.h'
       include 'pthpar.h'
       save
       integer ipath, i, j, k, inpath, nqmax, jfeff
       integer u2ipth, ntmp, nkpar, icode(micode)
       double precision reff, degen, xtmp, getsca
       double precision chi_r(maxpts),chi_i(maxpts)
       double precision tpakar(maxpts),tpaamp(maxpts),tpapha(maxpts)
       external         u2ipth, getsca

       reff      = -one
       xafs_path = 0
       inpath    = u2ipth(ipath)
cc       print*, ' xafs path ', ipath, u2ipth(ipath)
       if (inpath.le.0) return
       if (jpthff(inpath).le.0) then
          write(tmpstr,'(1x,a,i5)') ' no FEFF file for path ', ipath
          call warn(2,tmpstr)
          return
       else
          jfeff = jpthff(inpath)
          reff  = refpth(jfeff)
          degen = degpth(jfeff)
          xtmp  = ipath
          call setsca('path_index', xtmp)
          call setsca('degen', degen)
          call setsca('reff',  reff)
       endif
       call synvar
c
c path OK to use
       xafs_path = 1
       nkpar     = 0

       do 340 i = 1, maxpts
          tpakar(i) = (i-1)* qgrid
          tpaamp(i) = one
          tpapha(i) = zero
 340   continue

c loop over path parameters to evaluate
       do 500 i = 1, mpthpr
c
c set default path params here
          tmparr(1) = zero
          if (i.eq.jfps02)  then
             tmparr(1) = one
          elseif (i.eq.jfpdeg) then
             tmparr(1) = degpth(jfeff)
          endif
          ntmp = 0
          do 420 k = 1, micode
             icode(k)  = icdpar(k,i,inpath)
 420      continue
          if (iprint.ge.12)
     $          call rpndmp(icode)
c  evaluate parameter if it was defined
          if (icode(1).ne.0) then
             call decod(icode, micode,
     $            consts, scalar, array, narray, nparr,
     $            maxsize_array, maxarr,  ntmp, tmparr)
          end if
c
c set the param value from tmparr:
          param(i) = tmparr(1)
          if (i.eq.jfpkar) then
             if (nkpar.gt.0)  nkpar = min(nkpar, ntmp)
             if (nkpar.eq.0)  nkpar = ntmp
             do 470 j = 1, ntmp
                tpakar(j) = tmparr(j)
 470         continue
          elseif (i.eq.jfpaar) then
             if (nkpar.gt.0)  nkpar = min(nkpar, ntmp)
             if (nkpar.eq.0)  nkpar = ntmp
             do 480 j = 1, ntmp
                tpaamp(j) = tmparr(j)
 480         continue
          elseif (i.eq.jfppar) then
             if (nkpar.gt.0)  nkpar = min(nkpar, ntmp)
             if (nkpar.eq.0)  nkpar = ntmp
             do 490 j = 1, ntmp
                tpapha(j) = tmparr(j)
 490         continue
          endif
 500   continue

c
c   get chi(k) for this path from feff and path parameters
c
       if ( (inpath.gt.0).and.(jfeff.gt.0))  then
          if (iprint.ge.9) then
             call echo('calling chipth:')
c      c             print*, '   jfeff, maxpts   =', jfeff, maxpts, mffpts
c      c             print*, '   reff, degen     =', reff, degen
c      c          print*, 'sum_paths: jfeff, nffpts   =',
c      c     $         jfeff,nffpts(jfeff)
          end if

cc       print*, 'Phase array: ', nkpar, tpapha(1), tpapha(3)

          call chipth(theamp(1,jfeff), thepha(1,jfeff),
     $         qfeff(1,jfeff), xlamb(1,jfeff), realp(1,jfeff),
     $         nffpts(jfeff), reff, nkpar, tpakar, tpaamp, tpapha,
     $         maxpts, chi_r, chi_i)

cc          print*, ' path ', ipath, inpath, param(jfpdeg)
c          print*, 'degen,s02,e0,ei =',degen,
c     $         param(jfps02), param(jfpe0), param(jfpei)
c          print*, 'dphas,delr,sig2 =',
c     $         param(jfppha),param(jfpdr),param(jfpss2)
c          print*, 'theamp(28,jfeff)=', theamp(28,jfeff)
c          print*, 'qfeff(28,jfeff) =',  qfeff(28,jfeff)

cc          print*, chi_r(1), chi_i(1), chi_r(3), chi_i(3)

       end if
       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:
locate instance a at "a.fpp"
locate constant b at "b.fpp"
locate instance c at "c.fpp"
locate constant d at "d.fpp"

topology T {

  connections C {

    a.out[b] -> c.in[d]

  }

}

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
module rschroed_other

! Other Schroedinger integrators, not directly used by dftatom, but available
! for reuse. This module contains various rk4 integrators and Adams
! predictor-corrector integrators (both for outward and inward integration).

use types, only: dp
use ode1d_other, only: rk4_integrate_inward2
use utils, only: stop_error
use ode1d, only: adams_extrapolation_outward, adams_interp_outward, &
    rk4_integrate
use ode1d_other, only: adams_extrapolation_inward, adams_interp_inward

implicit none
private
public schroed_inward_adams_pr, schroed_outward_adams_pr, &
    integrate_rschroed_u_inward_rk4, integrate_rschroed_u_rk4

contains

subroutine integrate_rschroed_u_rk4(l, E, R, V, Vmid, P, Q, imax)
! Integrates the Schrodinger equation outward using Runge-Kutta 4th order method
!
! It integrates the Schroedinger equation in the P(r), Q(r) form outwards:
! P' = Q
! Q' = P'' = (2*(V-E) + l*(l+1)/R**2) * P
!
! Returns P(r), Q(r).
integer, intent(in) :: l
real(dp), intent(in) :: E
real(dp), intent(in) :: R(:)
real(dp), intent(in) :: V(:), Vmid(:)
real(dp), intent(out) :: P(:), Q(:)
integer, intent(out) :: imax
real(dp), parameter :: max_val = 1e300_dp

real(dp) :: y(2)
integer :: nr

real(dp), dimension(size(R)) :: C1, C2
real(dp), dimension(size(R)-1) :: C1mid, C2mid
real(dp), dimension(size(R)-1) :: Rmid

! y(:) are values of the components (2) of the equation, in our case:
! y(1) = P, y(2) = Q
! dydx(:) are the derivatives of the components (2) of the equation

nr = size(R)

y(1) = R(1)**(l+1)
y(2) = (l+1)*R(1)**l

C1 = 2*(V-E) + l*(l+1)/R**2
C2 = 0
Rmid = (R(:size(R)-1) + R(2:)) / 2
C1mid = 2*(Vmid-E) + l*(l+1)/Rmid**2
C2mid = 0
call rk4_integrate(R, y, C1, C2, C1mid, C2mid, max_val, P, Q, imax)
end subroutine

subroutine integrate_rschroed_u_inward_rk4(n, l, E, R, V, Vmid, &
        P, Q, imin)
! Integrates the Schrodinger equation inward using Runge-Kutta 4th order method
!
! It integrates the Schroedinger equation in the P(r), Q(r) form inwards:
! P' = Q
! Q' = P'' = (2*(V-E) + l*(l+1)/R**2) * P
!
! Returns P(r), Q(r).
integer, intent(in) :: l, n
real(dp), intent(in) :: E
real(dp), intent(in) :: R(:)
real(dp), intent(in) :: V(:), Vmid(:)
real(dp), intent(out) :: P(:), Q(:)
integer, intent(out) :: imin

real(dp) :: y(2), lambda, zeta, sigma
integer :: nr
integer :: i, i_max
real(dp), parameter :: max_val = 1e300_dp
integer, parameter :: asymptotic = 4
real(dp), dimension(size(R)) :: C1
real(dp), dimension(size(R)-1) :: C1mid
real(dp), dimension(size(R)-1) :: Rmid

! y(:) are values of the components (2) of the equation, in our case:
! y(1) = P, y(2) = Q
! dydx(:) are the derivatives of the components (2) of the equation

nr = size(R)

if (E > 0) then
    print *, "E =", E
    call stop_error("E < 0 required")
end if
if (nr < 4) call stop_error("integrate_rschroed_u_inward: n >= 4 required")
i_max = nr
lambda = sqrt(-2*E)
! One possible condition:
!do while (lambda*R(i_max) > 40)
! Start integrating at the 100x classical turning point
do while (R(i_max) > 100 * R(1))
    if (i_max == 2) then
        print *, "E =", E, "lambda =", lambda
        call stop_error("Can't start the inward integration")
    end if
    i_max = i_max - 1
end do
! The simplest asymptotic seems to be working just fine, but just in case, one
! can try more precise ones by changing the 'asymptotic' parameter above
if (asymptotic == 1) then
    ! Simplest asymptotic (simplified derivative)
    y(1) = R(i_max)**n * exp(-lambda*R(i_max))
    y(2) = - lambda * R(i_max)**n * exp(-lambda*R(i_max))
else if (asymptotic == 2) then
    ! Simple asymptotic (correct derivative)
    y(1) = R(i_max)**n * exp(-lambda*R(i_max))
    y(2) = (n * R(i_max)**(n-1) - lambda * R(i_max)**n) * exp(-lambda*R(i_max))
else if (asymptotic == 3) then
    ! More precise asymptotic:
    zeta = -V(i_max)*R(i_max)
    sigma = zeta/lambda
    y(1) = 0
    y(2) = 0
    ! Use 4 terms (we could use more if needed):
    do i = 0, 4
        y(1) = y(1) + R(i_max)**(sigma-i) * exp(-lambda*R(i_max)) * &
            a_k(i, lambda, sigma, l)
        y(2) = y(2) + R(i_max)**(sigma-i) * exp(-lambda*R(i_max)) * &
            b_k(i, lambda, sigma, l)
    end do
else if (asymptotic == 4) then
    lambda = sqrt(l*(l+1)/R(i_max)**2 + 2 * (V(i_max) - E))
    y(1) = 1
    y(2) = - lambda * R(i_max) * y(1)
else
    call stop_error("Unknown 'asymptotic' parameter")
end if

P(i_max+1:) = 0
Q(i_max+1:) = 0

C1 = 2*(V-E) + l*(l+1)/R**2
Rmid = (R(:size(R)-1) + R(2:)) / 2
C1mid = 2*(Vmid-E) + l*(l+1)/Rmid**2
call rk4_integrate_inward2(R(:i_max), y, C1(:i_max), C1mid(:i_max-1), &
        max_val, P(:i_max), Q(:i_max), imin)
end subroutine

real(dp) recursive function a_k(k, lambda, sigma, l) result(r)
integer, intent(in) :: k, l
real(dp), intent(in) :: lambda, sigma
if (k == 0) then
    r = 1
else
    r = (l*(l+1) - (sigma-k)*(sigma-k+1))/(2*k*lambda)
    r = r * a_k(k-1, lambda, sigma, l)
end if
end function

real(dp) recursive function b_k(k, lambda, sigma, l) result(r)
integer, intent(in) :: k, l
real(dp), intent(in) :: lambda, sigma
if (k == 0) then
    r = -lambda
else
    r = (-l*(l+1) + (sigma+k)*(sigma-k+1))/(2*k)
    r = r * a_k(k-1, lambda, sigma, l)
end if
end function

subroutine schroed_outward_adams_pr(l, E, R, Rp, V, P, Q, imax)
! Integrates the Schrodinger equation outward using Adams 4th order method
!
! It integrates the Schroedinger equation in the P(r), Q(r) form outwards using
! predictor-corrector:
! P' = Q
! Q' = P'' = C * P
! where C = 2*(V-E) + l*(l+1)/R**2
!
! Returns P(r), Q(r).
!
! It tranforms the problem to a uniform mesh 1 <= t <= N + 1, defined by
! R(t) and Rp(t) = dR/dt and Rpp(t) = d^2R/dt^2:
!
! u(t)   = P(R(t))
! up(t)  = du/dt  = dP/dR * dR/dt             = Q * Rp
! upp(t) = dup/dt = dQ/dR * Rp^2 + Q * dRp/dt = C*P*Rp^2 + Q * Rpp
!
! So
!
! up  = u * Rp
! upp = u * C * Rp^2 + up * Rpp / Rp
!
! For example if Rp = al*R, Rpp = al^2 * R, we get:
!
! up  = u * al * R
! upp = u * C * al^2 * R^2 + up * al

integer, intent(in) :: l
real(dp), intent(in) :: E
real(dp), intent(in) :: R(:), Rp(:)
real(dp), intent(in) :: V(:)
real(dp), intent(out) :: P(:), Q(:)
real(dp), parameter :: max_val = 1e6_dp
integer, intent(out) :: imax ! The integration was carried to R(imax)

real(dp), dimension(size(R)) :: C, u1, u2, u1p, u2p
integer :: i, it

C = (l*(l+1)/R**2 + 2 * (V-E))
u1(1:4) = R(1:4) ** (l+1)
u2(1:4) = (l+1) * R(1:4) ** l
u1p(1:4) = Rp(1:4)          * u2(1:4)
u2p(1:4) = Rp(1:4) * C(1:4) * u1(1:4)

do i = 4, size(R)-1
    u1(i+1)  = u1(i)  + adams_extrapolation_outward(u1p, i)
    u2(i+1)  = u2(i)  + adams_extrapolation_outward(u2p, i)
    do it = 1, 2
        u1p(i+1) = Rp(i+1)          * u2(i+1)
        u2p(i+1) = Rp(i+1) * C(i+1) * u1(i+1)
        u1(i+1)  = u1(i) + adams_interp_outward(u1p, i)
        u2(i+1)  = u2(i) + adams_interp_outward(u2p, i)
    end do
    if (abs(u1(i+1)) >= max_val .or. abs(u2(i+1)) >= max_val) then
        P = u1
        Q = u2
        imax = i
        return
    end if
end do
P = u1
Q = u2
imax = size(R)
end subroutine

subroutine schroed_inward_adams_pr(l, E, R, Rp, V, P, Q, imin)
! Integrates the Schrodinger equation inward using Adams 4th order method
!
! It integrates the Schroedinger equation in the P(r), Q(r) form inwards using
! predictor-corrector:
! P' = Q
! Q' = P'' = (2*(V-E) + l*(l+1)/R**2) * P
!
! Returns P(r), Q(r).
!
integer, intent(in) :: l
real(dp), intent(in) :: E
real(dp), intent(in) :: R(:), Rp(:)
real(dp), intent(in) :: V(:)
real(dp), intent(out) :: P(:), Q(:)
real(dp), parameter :: max_val = 1e300_dp
integer, intent(out) :: imin

real(dp), dimension(size(R)) :: C, u1, u2, u1p, u2p
real(dp) :: xkap
integer :: i, it, i_max

C = (l*(l+1)/R**2 + 2 * (V-E))

i_max = size(R)-4
do while (R(i_max) > 100 * R(1))
    if (i_max == 2) then
        call stop_error("Can't start the inward integration")
    end if
    i_max = i_max - 1
end do

xkap = sqrt(l*(l+1) / r(i_max)**2 + 2 * (V(i_max) - E))
u1(i_max:i_max+4) = exp(-xkap * (R(i_max:i_max+4) - R(i_max)))
u2(i_max:i_max+4) = - xkap * u1(i_max:i_max+4)
u1p(i_max:i_max+4) = Rp(i_max:i_max+4)                * u2(i_max:i_max+4)
u2p(i_max:i_max+4) = Rp(i_max:i_max+4) * C(i_max:i_max+4) * u1(i_max:i_max+4)

do i = i_max, 2, -1
    u1(i-1)  = u1(i)  + adams_extrapolation_inward(u1p, i)
    u2(i-1)  = u2(i)  + adams_extrapolation_inward(u2p, i)
    do it = 1, 2
        u1p(i-1) = Rp(i-1)          * u2(i-1)
        u2p(i-1) = Rp(i-1) * C(i-1) * u1(i-1)
        u1(i-1)  = u1(i) + adams_interp_inward(u1p, i)
        u2(i-1)  = u2(i) + adams_interp_inward(u2p, i)
    end do
    if (abs(u1(i-1)) >= max_val .or. abs(u2(i-1)) >= max_val) then
        P = u1
        Q = u2
        P(i_max+1:) = 0
        Q(i_max+1:) = 0
        imin = i
        return
    end if
end do
P = u1
Q = u2
P(i_max+1:) = 0
Q(i_max+1:) = 0
imin = 1
end subroutine

end module

max_stars_repo_name: rc/dftatom

max_issues_repo_name: rc/dftatom

max_forks_repo_name: rc/dftatom


!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.

!> @brief Reads vs30 and thickness from VS_RS.out
!! @author Ilario Mazzieri
!> @date September, 2015
!> @version 1.0
!> @param[in] filec file name for the mesh
!> @param[in] num_tria  number of triangular elements
!> @param[out] velocity vs30
!> @param[out] thickness thickness of soft sediments

      subroutine READ_FILEVS(filec,num_tria,velocity,sediments)

      implicit none

      integer*4 :: num_tria
      integer*4 :: i
      integer*4 :: trash

      character*70 :: filec
      real*8, dimension(num_tria) :: velocity, sediments

      open(20,file=filec)

      do i = 1,num_tria
        read(20,*) trash, velocity(i), sediments(i)
      enddo

      close(20)

      return

      end subroutine READ_FILEVS


text:
      program opkdemo5
c-----------------------------------------------------------------------
c Demonstration program for DLSODPK.
c ODE system from ns-species interaction pde in 2 dimensions.
c This is the version of 14 June 2001.
c
c This version is in double precision.
c-----------------------------------------------------------------------
c This program solves a stiff ODE system that arises from a system
c of partial differential equations.  The PDE system is a food web
c population model, with predator-prey interaction and diffusion on
c the unit square in two dimensions.  The dependent variable vector is
c
c         1   2        ns
c   c = (c , c , ..., c  )
c
c and the PDEs are as follows:
c
c     i               i      i
c   dc /dt  =  d(i)*(c    + c   )  +  f (x,y,c)  (i=1,...,ns)
c                     xx     yy        i
c
c where
c                  i          ns         j
c   f (x,y,c)  =  c *(b(i) + sum a(i,j)*c )
c    i                       j=1
c
c The number of species is ns = 2*np, with the first np being prey and
c the last np being predators.  The coefficients a(i,j), b(i), d(i) are:
c
c   a(i,i) = -a  (all i)
c   a(i,j) = -g  (i .le. np, j .gt. np)
c   a(i,j) =  e  (i .gt. np, j .le. np)
c   b(i) =  b*(1 + alpha*x*y)  (i .le. np)
c   b(i) = -b*(1 + alpha*x*y)  (i .gt. np)
c   d(i) = dprey  (i .le. np)
c   d(i) = dpred  (i .gt. np)
c
c The various scalar parameters are set in subroutine setpar.
c
c The boundary conditions are: normal derivative = 0.
c A polynomial in x and y is used to set the initial conditions.
c
c The PDEs are discretized by central differencing on a mx by my mesh.
c
c The ODE system is solved by DLSODPK using method flag values
c mf = 10, 21, 22, 23, 24, 29.  The final time is tmax = 10, except
c that for mf = 10 it is tmax = 1.0d-3 because the problem is stiff,
c and for mf = 23 and 24 it is tmax = 2 because the lack of symmetry
c in the problem makes these methods more costly.
c
c Two preconditioner matrices are used.  One uses a fixed number of
c Gauss-Seidel iterations based on the diffusion terms only.
c The other preconditioner is a block-diagonal matrix based on
c the partial derivatives of the interaction terms f only, using
c block-grouping (computing only a subset of the ns by ns blocks).
c For mf = 21 and 22, these two preconditioners are applied on
c the left and right, respectively, and for mf = 23 and 24 the product
c of the two is used as the one preconditioner matrix.
c For mf = 29, the inverse of the product is applied.
c
c Two output files are written: one with the problem description and
c and performance statistics on unit 6, and one with solution profiles
c at selected output times (for mf = 22 only) on unit 8.
c-----------------------------------------------------------------------
c Note: In addition to the main program and 10 subroutines
c given below, this program requires the LINPACK subroutines
c DGEFA and DGESL, and the BLAS routine DAXPY.
c-----------------------------------------------------------------------
c Reference:
c     Peter N. Brown and Alan C. Hindmarsh,
c     Reduced Storage Matrix Methods in Stiff ODE Systems,
c     J. Appl. Math. & Comp., 31 (1989), pp. 40-91;
c     Also LLNL Report UCRL-95088, Rev. 1, June 1987.
c-----------------------------------------------------------------------
      external fweb, jacbg, solsbg
      integer ns, mx, my, mxns,
     1        mp, mq, mpsq, itmax,
     2        meshx,meshy,ngx,ngy,ngrp,mxmp,jgx,jgy,jigx,jigy,jxr,jyr
      integer i, imod3, iopt, iout, istate, itask, itol, iwork,
     1   jacflg, jpre, leniw, lenrw, liw, lrw, mf,
     2   ncfl, ncfn, neq, nfe, nfldif, nfndif, nli, nlidif, nni, nnidif,
     3   nout, npe, nps, nqu, nsdif, nst
      double precision aa, ee, gg, bb, dprey, dpred,
     1     ax, ay, acoef, bcoef, dx, dy, alph, diff, cox, coy,
     2     uround, srur
      double precision avdim, atol, cc, hu, rcfl, rcfn, dumach,
     1   rtol, rwork, t, tout
      dimension neq(1), atol(1), rtol(1)
c
c The problem Common blocks below allow for up to 20 species,
c up to a 50x50 mesh, and up to a 20x20 group structure.
      common /pcom0/ aa, ee, gg, bb, dprey, dpred
      common /pcom1/ ax, ay, acoef(20,20), bcoef(20), dx, dy, alph,
     1               diff(20), cox(20), coy(20), ns, mx, my, mxns
      common /pcom2/ uround, srur, mp, mq, mpsq, itmax
      common /pcom3/ meshx, meshy, ngx, ngy, ngrp, mxmp,
     2   jgx(21), jgy(21), jigx(50), jigy(50), jxr(20), jyr(20)
c
c The dimension of cc below must be .ge. 2*neq, where neq = ns*mx*my.
c The dimension lrw of rwork must be .ge. 17*neq + ns*ns*ngrp + 61,
c and the dimension liw of iwork must be .ge. 35 + ns*ngrp.
      dimension cc(576), rwork(5213), iwork(67)
      data lrw/5213/, liw/67/
c
      open (unit=6, file='demout', status='new')
      open (unit=8, file='ccout', status='new')
c
      ax = 1.0d0
      ay = 1.0d0
c
c Call setpar to set problem parameters.
      call setpar
c
c Set remaining problem parameters.
      neq(1) = ns*mx*my
      mxns = mx*ns
      dx = ax/(mx-1)
      dy = ay/(my-1)
      do 10 i = 1,ns
        cox(i) = diff(i)/dx**2
 10     coy(i) = diff(i)/dy**2
c
c Write heading.
      write(6,20)ns, mx,my,neq(1)
 20   format(' Demonstration program for DLSODPK package'//
     1   ' Food web problem with ns species, ns =',i4/
     2   ' Predator-prey interaction and diffusion on a 2-d square'//
     3   ' Mesh dimensions (mx,my) =',2i4/
     4   ' Total system size is neq =',i7//)
      write(6,25) aa,ee,gg,bb,dprey,dpred,alph
 25   format(' Matrix parameters:  a =',d12.4,'   e =',d12.4,
     1   '   g =',d12.4/20x,' b =',d12.4//
     2   ' Diffusion coefficients: dprey =',d12.4,'   dpred =',d12.4/
     3   ' Rate parameter alpha =',d12.4//)
c
c Set remaining method parameters.
      jpre = 3
      jacflg = 1
      iwork(3) = jpre
      iwork(4) = jacflg
      iopt = 0
      mp = ns
      mq = mx*my
      mpsq = ns*ns
      uround = dumach()
      srur = sqrt(uround)
      meshx = mx
      meshy = my
      mxmp = meshx*mp
      ngx = 2
      ngy = 2
      ngrp = ngx*ngy
      call gset (meshx, ngx, jgx, jigx, jxr)
      call gset (meshy, ngy, jgy, jigy, jyr)
      iwork(1) = mpsq*ngrp
      iwork(2) = mp*ngrp
      itmax = 5
      itol = 1
      rtol(1) = 1.0d-5
      atol(1) = rtol(1)
      itask = 1
      write(6,30)ngrp,ngx,ngy,itmax,rtol(1),atol(1)
 30   format(' Preconditioning uses interaction-only block-diagonal',
     1   ' matrix'/' with block-grouping, and Gauss-Seidel iterations'//
     2   ' Number of diagonal block groups = ngrp =',i4,
     3   '   (ngx by ngy, ngx =',i2,'  ngy =',i2,' )'//
     4   ' G-S preconditioner uses itmax iterations, itmax =',i3//
     5   ' Tolerance parameters: rtol =',d10.2,'   atol =',d10.2)
c
c
c Loop over mf values 10, 21, ..., 24, 29.
c
      do 90 mf = 10,29
      if (mf .gt. 10 .and. mf .lt. 21) go to 90
      if (mf .gt. 24 .and. mf .lt. 29) go to 90
      write(6,40)mf
 40   format(//80('-')//' Solution with mf =',i3//
     1   '   t       nstep  nfe  nni  nli  npe  nq',
     2   4x,'h          avdim    ncf rate    lcf rate')
c
      t = 0.0d0
      tout = 1.0d-8
      nout = 18
      if (mf .eq. 10) nout = 6
      if (mf .eq. 23 .or. mf .eq. 24) nout = 10
      call cinit (cc)
      if (mf .eq. 22) call outweb (t, cc, ns, mx, my, 8)
      istate = 1
      nli = 0
      nni = 0
      ncfn = 0
      ncfl = 0
      nst = 0
c
c Loop over output times and call DLSODPK.
c
      do 70 iout = 1,nout
        call dlsodpk (fweb, neq, cc, t, tout, itol, rtol, atol, itask,
     1         istate, iopt, rwork, lrw, iwork, liw, jacbg, solsbg, mf)
        nsdif = iwork(11) - nst
        nst = iwork(11)
        nfe = iwork(12)
        npe = iwork(13)
        nnidif = iwork(19) - nni
        nni = iwork(19)
        nlidif = iwork(20) - nli
        nli = iwork(20)
        nfndif = iwork(22) - ncfn
        ncfn = iwork(22)
        nfldif = iwork(23) - ncfl
        ncfl = iwork(23)
        nqu = iwork(14)
        hu = rwork(11)
        avdim = 0.0d0
        rcfn = 0.0d0
        rcfl = 0.0d0
        if (nnidif .gt. 0) avdim = real(nlidif)/real(nnidif)
        if (nsdif .gt. 0) rcfn = real(nfndif)/real(nsdif)
        if (nnidif .gt. 0) rcfl = real(nfldif)/real(nnidif)
        write(6,50)t,nst,nfe,nni,nli,npe,nqu,hu,avdim,rcfn,rcfl
 50     format(d10.2,i5,i6,3i5,i4,2d11.2,d10.2,d12.2)
        imod3 = iout - 3*(iout/3)
        if (mf .eq. 22 .and. imod3 .eq. 0) call outweb (t,cc,ns,mx,my,8)
        if (istate .eq. 2) go to 65
        write(6,60)t
 60     format(//' final time reached =',d12.4//)
        go to 75
 65     continue
        if (tout .gt. 0.9d0) tout = tout + 1.0d0
        if (tout .lt. 0.9d0) tout = tout*10.0d0
 70     continue
c
 75   continue
      nst = iwork(11)
      nfe = iwork(12)
      npe = iwork(13)
      lenrw = iwork(17)
      leniw = iwork(18)
      nni = iwork(19)
      nli = iwork(20)
      nps = iwork(21)
      if (nni .gt. 0) avdim = real(nli)/real(nni)
      ncfn = iwork(22)
      ncfl = iwork(23)
      write (6,80) lenrw,leniw,nst,nfe,npe,nps,nni,nli,avdim,
     1               ncfn,ncfl
 80   format(//' Final statistics for this run:'/
     1   ' rwork size =',i8,'   iwork size =',i6/
     2   ' number of time steps            =',i5/
     3   ' number of f evaluations         =',i5/
     4   ' number of preconditioner evals. =',i5/
     4   ' number of preconditioner solves =',i5/
     5   ' number of nonlinear iterations  =',i5/
     5   ' number of linear iterations     =',i5/
     6   ' average subspace dimension  =',f8.4/
     7   i5,' nonlinear conv. failures,',i5,' linear conv. failures')
c
 90   continue
c      stop
c------  end of main program for DLSODPK demonstration program ----------
      end

ext: f

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
      subroutine gen9_rap(r,p,wt9,*)
      implicit none
      include 'constants.f'
      include 'masses.f'
      include 'mxdim.f'
      include 'breit.f'
      integer nu
      double precision r(mxdim)
      double precision wt9,p(mxpart,4),tp(4),tm(4),bp(4),bm(4),
     & n(4),e(4),g(4),ep(4),em(4),nn(4),nb(4),wp(4),wm(4),
     & wtepnn,wtnbem
      double precision wtttwg,wtwp,wtwm,s3min,wt0
      parameter(wt0=1d0/twopi**4)

*     q(-p1) +qbar(-p2)=t(nu(p3)+e^+(p4)+b(p5))                        *
*                       +t~(b~(p6)+e^-(p7)+nu(p8))+g(p9)               *
*                                                                      *

      wt9=0d0
C---call gen5 that uses r(1)....r(13)
      call gen5(r,p,wtttwg,*999)
      wtttwg=(mt*twidth*pi)**2*wtttwg
      do nu=1,4
      n(nu)=p(3,nu)
      e(nu)=p(4,nu)
      tp(nu)=p(5,nu)
      tm(nu)=p(6,nu)
      g(nu)=p(7,nu)
      enddo

      s3min=0d0
      n3=1
      mass3=wmass
      width3=wwidth

      call phi1_2m(mb,r(14),r(15),r(16),s3min,tp,bm,wp,wtwp,*999)
      call phi1_2m(mb,r(17),r(18),r(19),s3min,tm,bp,wm,wtwm,*999)
      call phi3m0(r(20),r(21),wp,nn,ep,wtepnn,*999)
      call phi3m0(r(22),r(23),wm,em,nb,wtnbem,*999)
      wt9=wt0*wtepnn*wtnbem*wtwp*wtwm*wtttwg
*     q(-p1) +qbar(-p2)=t(nu(p3)+e^+(p4)+b(p5))                        *
*                       +t~(b~(p6)+e^-(p7)+nu(p8)+W(nu(p9)+e+(p10)     *
*                                                                      *
      do nu=1,4
      p(3,nu)=nn(nu)
      p(4,nu)=ep(nu)
      p(5,nu)=bm(nu)
      p(6,nu)=bp(nu)
      p(7,nu)=em(nu)
      p(8,nu)=nb(nu)
      p(9,nu)=n(nu)
      p(10,nu)=e(nu)
      p(11,nu)=g(nu)
      enddo

      return
 999  wt9=0d0
      p(:,:)=0d0
      return 1

      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
c     UC Copyright Notice
c
c     This software is Copyright (c) 2014-2015 The Regents of the
c     University of California. All Rights Reserved.
c
c     Permission to copy and modify this software and its documentation
c     for educational, research and non-profit purposes, without fee,
c     and without a written agreement is hereby granted, provided that
c     the above copyright notice, this paragraph and the following three
c     paragraphs appear in all copies.
c
c     Permission to make commercial use of this software may be obtained
c     by contacting:
c
c     Technology Transfer Office
c     9500 Gilman Drive, Mail Code 0910
c     University of California
c     La Jolla, CA 92093-0910
c     (858) 534-5815
c     invent@ucsd.edu
c
c     This software program and documentation are copyrighted by The
c     Regents of the University of California. The software program and
c     documentation are supplied "as is", without any accompanying
c     services from The Regents. The Regents does not warrant that the
c     operation of the program will be uninterrupted or error-free. The
c     end-user understands that the program was developed for research
c     purposes and is advised not to rely exclusively on the program for
c     any reason.
c
c     IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY
c     PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL
c     DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS
c     SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF
c     CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
c     THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY
c     WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
c     OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE
c     SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE
c     UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
c     MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.

c Portions of the code Copyright (c) 2009-2011 Open Source Medical
c Software Corporation, University of California, San Diego.
c
c Portions of the code Copyright (c) 2000-2007, Stanford University,
c     Rensselaer Polytechnic Institute, Kenneth E. Jansen,
c     Charles A. Taylor.
c
c  See SimVascular Acknowledgements file for additional
c  contributors to the source code.
c
c  Redistribution and use in source and binary forms, with or without
c  modification, are permitted provided that the following conditions
c  are met:
c
c  Redistributions of source code must retain the above copyright notice,
c  this list of conditions and the following disclaimer.
c  Redistributions in binary form must reproduce the above copyright
c  notice, this list of conditions and the following disclaimer in the
c  documentation and/or other materials provided with the distribution.
c  Neither the name of the Stanford University or Rensselaer Polytechnic
c  Institute nor the names of its contributors may be used to endorse or
c  promote products derived from this software without specific prior
c  written permission.
c
c  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
c  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
c  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
c  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
c  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
c  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
c  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
c  OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
c  AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
c  OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
c  THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
c  DAMAGE.
c
c

#include "cvFlowsolverOptions.h"

!> This routine computes and assembles the data corresponding to the
!! boundary elements.

      subroutine AsBNASC ( x,       shpb,
     &                   ienb,  iBCB)
        use pvsQbi
        include "global.h"
        include "common_blocks/conpar.h"
        include "common_blocks/elmpar.h"
        include "common_blocks/genpar.h"
        include "common_blocks/intpt.h"
        include "common_blocks/propar.h"
        include "common_blocks/shpdat.h"
        include "common_blocks/nomodule.h"
c
c..............................Declaration...............................................
       INTEGER ienb,iBCB,iel,iface,ipt3,k,lnode,n,nodlcl
       REAL*8 x,shpb,bnorm,dxdxib,rl,sgn,shdrv,shglb,shpfun,temp,temp1,
     &temp2,temp3,v1,v2
       REAL*8 WdetJb,xlb
c........................................................................................

        dimension xlb(npro,nenl,nsd),
     &            rl(npro,nshl),     WdetJb(npro)

        dimension x(numnp,nsd),
     &            shpb(nshl,ngaussb),      shglb(nsd,nshl,ngaussb),
     &            ienb(npro,nshl),
     &            iBCB(npro,ndiBCB)

        dimension lnode(27),               sgn(npro,nshl),
     &            shpfun(npro,nshl),        shdrv(npro,nsd,nshl)

c
        dimension dxdxib(npro,nsd,nsd),      temp(npro),
     &            temp1(npro),               temp2(npro),
     &            temp3(npro),
     &            v1(npro,nsd),              v2(npro,nsd)

#if (VER_CLOSEDLOOP == 1)
      LOGICAL flagDirichlet
      INTEGER i
#endif

c
c.... get the matrix of mode signs for the hierarchic basis functions
c
        if (ipord .gt. 1) then
           call getsgn(ienb,sgn)
        endif
c
c.... gather the variables
c
        call localx(x,      xlb,    ienb,   nsd,    'gather  ')
c
c.... get the boundary element residuals
c
        rl  = zero
c
c.... compute the nodes which lie on the boundary (hierarchic)
c
        call getbnodes(lnode)
c
c.... loop through the integration points
c
        ngaussb = nintb(lcsyst)

        do intp = 1, ngaussb
c
c.... get the hierarchic shape functions at this int point
c
           shglb=zero  ! protect debugger
           call getshpb(shpb,        shglb,        sgn,
     &              shpfun,       shdrv)

c
c.... compute the normal to the boundary. This is achieved by taking
c     the cross product of two vectors in the plane of the 2-d
c     boundary face.
c
           if(lcsyst.ne.6) then
              ipt3=3
           else
              ipt3=5
           endif
           v1 = xlb(:,2,:) - xlb(:,1,:)
           v2 = xlb(:,ipt3,:) - xlb(:,1,:)

c
c.....The following are done in order to correct temp1..3
c     based on the results from compressible code.  This is done only
c     for wedges, depending on the boundary face.(tri or quad)
c
           if (lcsyst .eq. 1) then
              temp1 = v1(:,2) * v2(:,3) - v2(:,2) * v1(:,3)
              temp2 = v2(:,1) * v1(:,3) - v1(:,1) * v2(:,3)
              temp3 = v1(:,1) * v2(:,2) - v2(:,1) * v1(:,2)
           else
              temp1 = - v1(:,2) * v2(:,3) + v2(:,2) * v1(:,3)
              temp2 = - v2(:,1) * v1(:,3) + v1(:,1) * v2(:,3)
              temp3 = - v1(:,1) * v2(:,2) + v2(:,1) * v1(:,2)
           endif
c
           temp       = one / sqrt ( temp1**2 + temp2**2 + temp3**2 )

           if (lcsyst .eq. 3) then
              WdetJb     = (1 - Qwtb(lcsyst,intp)) / (four*temp)
           elseif (lcsyst .eq. 4) then
              WdetJb     = Qwtb(lcsyst,intp) / temp
           else
              WdetJb     = Qwtb(lcsyst,intp) / (four*temp)
           endif

c......here I only want the d Gamma, not the n_i
           do iel=1,npro
              if (btest(iBCB(iel,1),1)) then
              else

#if (VER_CLOSEDLOOP == 1)
                flagDirichlet = .TRUE.
                DO i=1, numDirichletSrfs
                  IF (iBCB(iel,2) .EQ. nsrflistDirichlet(i)) THEN
                    flagDirichlet = .FALSE.
                  END IF
                END DO
                IF (flagDirichlet) THEN
                  WdetJb(iel) = zero  ! we want zeros where we are not integrating
                END IF
#else
                WdetJb(iel) = zero  ! we want zeros where we are not integrating
#endif

              endif
           enddo

c
c  Now lets calculate Integral N_(a:e) d Gamma
c
c
           do n = 1, nshlb
              nodlcl = lnode(n)
            rl(:,nodlcl) = rl(:,nodlcl) + shpfun(:,nodlcl) * WdetJb(:)

           enddo

        enddo  ! quadrature point loop
c
c.... assemble the NASC vector
c
        call local (NASC,    rl,     ienb,   1,  'scatter ')
c
c     push the surf number which we have associated with boundary
C     elements up to the global level in the array ndsurf
c
        do iel=1,npro
           if (iBCB(iel,2) .ne. 0) then
              iface = iBCB(iel,2)
              ndsurf(ienb(iel,1:nshlb))=iface
           endif
        enddo
c
c.... end
c
        return
        end


text:
C------------------------------------------------------------------------------
        subroutine ftgics(iunit,xrval,yrval,xrpix,yrpix,xinc,yinc,rot,
     &                   type,status)

C       read the values of the celestial coordinate system keywords.
C       These values may be used as input to the subroutines that
C       calculate celestial coordinates. (FTXYPX, FTWLDP)

C       This routine assumes that the CHDU contains an image
C       with the RA type coordinate running along the first axis
C       and the DEC type coordinate running along the 2nd axis.

        double precision xrval,yrval,xrpix,yrpix,xinc,yinc,rot
        integer iunit,status,tstat
        character*(*) type
        character comm*20,ctype*8

        if (status .gt. 0)return

        call ftgkyd(iunit,'CRVAL1',xrval,comm,status)
        call ftgkyd(iunit,'CRVAL2',yrval,comm,status)

        call ftgkyd(iunit,'CRPIX1',xrpix,comm,status)
        call ftgkyd(iunit,'CRPIX2',yrpix,comm,status)

        call ftgkyd(iunit,'CDELT1',xinc,comm,status)
        call ftgkyd(iunit,'CDELT2',yinc,comm,status)

        call ftgkys(iunit,'CTYPE1',ctype,comm,status)

        if (status .gt. 0)then
            call ftpmsg('FTGICS could not find all the required'//
     &                  'celestial coordinate Keywords.')
            status=505
            return
        end if

        type=ctype(5:8)

        tstat=status
        call ftgkyd(iunit,'CROTA2',rot,comm,status)
        if (status .gt. 0)then
C           CROTA2 is assumed to = 0 if keyword is not present
            status=tstat
            rot=0.
        end if
        end

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
!
! Module for normal score transform related routines
! Author: John Manchuk, 2010
!
! Module Dependencies:
!   NORMALDIST - used to compute inverse normal distribution
!   QUICKSORT  - for sorting transform tables
!   ACORN_RNG  - random number generator for breaking ties
!
module normalscore
    use normaldist, ONLY : snorm_inv, snorm_cdf, locate, powint
    use quicksort
    use acorn_rng
    implicit none

    !These are general functions.  All other functions in this module are
    !called from these depending on the inputs.
    public :: build_transtable, normal_score

    type :: transtable !used for jagged arrays
        integer :: n
        real*8 :: sumw
        real*8 :: mint
        real*8, allocatable :: table(:,:)
    end type transtable
    public :: transtable

    !Machine precision
    real*8, parameter :: MACHP = sqrt(epsilon(real(1,8)))

    !Random number generator for breaking ties
    integer, parameter :: seed = 69069
    logical :: isacorninit = .false.
    type(acorn_type) :: rng

    !Work array
    real*8, allocatable :: rwork(:,:)

    contains

    !Compute the normal scores
    subroutine normal_score ( data_, colvr, colct, ct, numq, tmin, tmax, ptt, nsvr, test )
        real*8, pointer :: data_(:,:), nsvr(:)
        integer, intent (in) :: colvr, colct !variable and category columns
        integer, intent (in) :: ct !the category of interest
        integer, intent (in) :: numq
        real*8, intent (in) :: tmin, tmax !trimming limits
        type(transtable), target, intent(in) :: ptt
        logical, intent(out) :: test
        integer :: ncol, nrow, nq
        integer :: i, j, maxiter
        real*8, pointer :: pvr(:),pvrns(:)

        !Initialize RNG
        if(.not.isacorninit)then
            call acorninit(rng,seed)
            isacorninit = .true.
        endif

        !Choose the version based on inputs
        if( colct < 1 ) then
            ncol = ubound(data_,1)
            nrow = ubound(data_,2)
            nq = ptt%n
            pvr => ptt%table(1,:)
            pvrns => ptt%table(2,:)
            do i = 1, nrow
                if(data_(colvr,i) < tmin .or. data_(colvr,i) >= tmax) then
                    nsvr(i) = -999.
                else
                    nsvr(i) = data_(colvr,i) + acorni(rng)*MACHP
                    maxiter = 0; test = .true.
                    do while ( nsvr(i) < pvr(1) .or. nsvr(i) > pvr(nq) )
                        nsvr(i) = data_(colvr,i) + acorni(rng)*MACHP
                        maxiter = maxiter + 1
                        if(maxiter > 1E7)then
                            test = .false.
                            return
                        endif
                    enddo
                    call locate(pvr,nq,1,nq,nsvr(i),j)
                    j = min(max(1,j),(nq-1))
                    nsvr(i) = powint(pvr(j),pvr(j+1),pvrns(j),pvrns(j+1),nsvr(i),1.D0)
                endif
            enddo
        else
            call normal_score_by_category ( data_, colvr, colct, ct, &
                numq, tmin, tmax, ptt, nsvr )
            test = .true.
        endif
        return
    end subroutine normal_score

    !Compute the normal scores of a category
    subroutine normal_score_by_category ( data_, colvr, colct, ct, numq, tmin, tmax, ptt, nsvr )
        real*8, pointer :: data_(:,:), nsvr(:)
        integer, intent (in) :: colvr, colct !variable and category columns
        integer, intent (in) :: ct !the category of interest
        real*8, intent (in) :: tmin, tmax !trimming limits
        type(transtable), target, intent(in) :: ptt
        integer :: ncol, nrow, numq
        integer :: i, j
        real*8, pointer :: pvr(:),pvrns(:)

        ncol = ubound(data_,1)
        nrow = ubound(data_,2)
        numq = ptt%n
        pvr => ptt%table(1,:)
        pvrns => ptt%table(2,:)

        do i = 1, nrow
            if(data_(colvr,i) < tmin .or. data_(colvr,i) >= tmax) then
                nsvr(i) = -999.
            else
                if(floor(data_(colct,i)+0.1) == ct)then
                    nsvr(i) = data_(colvr,i) + acorni(rng)*MACHP
                    call locate(pvr,numq,1,numq,nsvr(i),j)
                    j = min(max(1,j),(numq-1))
                    nsvr(i) = powint(pvr(j),pvr(j+1),pvrns(j),pvrns(j+1),nsvr(i),1.D0)
                endif
            endif
        enddo
        return
    end subroutine normal_score_by_category

    !build a transformation table
    logical function build_transtable ( data_, colvr, colwt, colct, ct, numq, tmin, tmax, ptt )
        real*8, pointer :: data_(:,:) !a block of data for the table, must be ncol by nrows
        integer, intent (in) :: colvr, colwt, colct !variable, weight, and category columns
        integer, intent (in) :: ct !the category of interest
        integer, intent (in) :: numq !number of quantiles
        real*8, intent (in) :: tmin, tmax !trimming limits
        type(transtable), intent(out) :: ptt !Transform table to be built
        integer :: ncol, nrow
        integer :: i, j, test
        integer :: numdata
        real*8 :: mindat

        !Some initial checks
        build_transtable = .false.
        ncol = ubound(data_,1)
        nrow = ubound(data_,2)
        if(colvr < 1)then
            write(*,'(a)') 'ERROR BUILDING TRANSTABLE, INVALID VARIABLE COLUMN'
            return
        endif
        if(colvr > ncol)then
            write(*,'(a)') 'ERROR BUILDING TRANSTABLE, REQUESTED VARIABLE COLUMN TOO LARGE'
            return
        endif
        if(colwt > ncol)then
            write(*,'(a)') 'ERROR BUILDING TRANSTABLE, REQUESTED WEIGHT COLUMN TOO LARGE'
            return
        endif
        if(colct > ncol)then
            write(*,'(a)') 'ERROR BUILDING TRANSTABLE, REQUESTED CATEGORY COLUMN TOO LARGE'
            return
        endif
        if(numq > 0 .and. numq < 3)then
            write(*,'(a)') 'ERROR BUILDING TRANSTABLE, NUMBER OF QUANTILES MUST BE > 2'
            return
        endif
        if(nrow < 2)then
            write(*,'(a)') 'ERROR BUILDING TRANSTABLE, NOT ENOUGH DATA'
            return
        endif
        if(allocated(ptt%table)) deallocate(ptt%table)

        call acorninit(rng,seed)
        numdata = 0

        !Choose the version of this routine depending on the inputs
        if(colwt < 1 .and. colct < 1)then
            !Count the data
            do i = 1, nrow
                if(data_(colvr,i) >= tmin .and. data_(colvr,i) < tmax)then
                    numdata = numdata + 1
                endif
            enddo
            if(numdata < 2) return

            !Fill in the transform table
            allocate(ptt%table(2,numdata), stat = test)
            if(test /= 0)then
                write(*,'(/,a,/)') 'ERROR BUILDING TRANSTABLE, ALLOCATION FAILED'
                return
            endif
            j = 0
            mindat = huge(mindat)
            do i = 1, nrow
                if(data_(colvr,i) >= tmin .and. data_(colvr,i) < tmax)then
                    j = j + 1
                    if(data_(colvr,i) < mindat) mindat = data_(colvr,i)
                    ptt%table(1,j) = data_(colvr,i) + acorni(rng)*MACHP
                    ptt%table(2,j) = 1.
                endif
            enddo
            ptt%mint = mindat

            !define size of the table
            ptt%n = numdata
            ptt%sumw = real(numdata,8)

            ! Compute nscores
            call nscore_transtable ( ptt, numq )

            build_transtable = .true.
        else
            if(colwt < 1)then !colct > 0
                build_transtable = build_transtable_by_category ( data_, &
                             colvr, colct, ct, numq, tmin, tmax, ptt )
            elseif(colct < 1)then !colwt > 0
                build_transtable = build_transtablew ( data_, &
                             colvr, colwt, numq, tmin, tmax, ptt )
            else !colwt and colct > 0
                build_transtable = build_transtablew_by_category ( data_, &
                             colvr, colwt, colct, ct, numq, tmin, tmax, ptt )
            endif
        endif

        return
    end function build_transtable

    !Build a weighted transform table
    ! this routine is called from build_transtable so error checking is not needed
    logical function build_transtablew ( data_, colvr, colwt, numq, tmin, tmax, ptt )
        real*8, pointer :: data_(:,:) !a block of data for the table, must be ncol by nrows
        integer, intent (in) :: colvr, colwt !variable and category columns
        integer, intent (in) :: numq !number of quantiles
        real*8, intent (in) :: tmin, tmax !trimming limits
        type(transtable), intent(inout) :: ptt !pointer to the transform table to be built
        integer :: ncol, nrow
        integer :: i, j, test
        integer :: numdata
        real*8 :: sumweight

        build_transtablew = .false.
        ncol = ubound(data_,1)
        nrow = ubound(data_,2)
        numdata = 0
        sumweight = 0.

        !Count the data
        do i = 1, nrow
            if(data_(colvr,i) >= tmin .and. data_(colvr,i) < tmax)then
                numdata = numdata + 1
            endif
        enddo
        if(numdata < 2) return

        !Fill in the transform table
        allocate(ptt%table(2,numdata), stat = test)
        if(test /= 0)then
            write(*,'(/,a,/)') 'ERROR BUILDING TRANSTABLE, ALLOCATION FAILED'
            return
        endif
        j = 0
        do i = 1, nrow
            if(data_(colvr,i) >= tmin .and. data_(colvr,i) < tmax )then
                j = j + 1
                ptt%table(1,j) = data_(colvr,i) + acorni(rng)*MACHP
                ptt%table(2,j) = data_(colwt,i)
                if(ptt%table(2,j) < MACHP )then ! Check if the weight is numerically small
                    ptt%table(2,j) = MACHP
                endif
                sumweight = sumweight + ptt%table(2,j)
            endif
        enddo

        !define size of the table
        ptt%n = numdata
        ptt%sumw = sumweight

        ! Compute nscores
        call nscore_transtable ( ptt, numq )

        build_transtablew = .true.
        return
    end function build_transtablew

    !Build a transform table by category
    ! this routine is called from build_transtable so error checking is not needed
    logical function build_transtable_by_category ( data_, colvr, colct, ct, numq, tmin, tmax, ptt )
        real*8, pointer :: data_(:,:) !a block of data for the table, must be ncol by nrows
        integer, intent (in) :: colvr, colct !variable and category columns
        integer, intent (in) :: ct !category of interest
        integer, intent (in) :: numq !number of quantiles
        real*8, intent (in) :: tmin, tmax !trimming limits
        type(transtable), intent(inout) :: ptt !pointer to the transform table to be built
        integer :: ncol, nrow
        integer :: i, j, test
        integer :: numdata

        build_transtable_by_category = .false.
        ncol = ubound(data_,1)
        nrow = ubound(data_,2)
        numdata = 0

        !Count the data
        do i = 1, nrow
            if(data_(colvr,i) >= tmin .and. data_(colvr,i) < tmax .and. &
               floor(data_(colct,i)+0.1) == ct)then
                numdata = numdata + 1
            endif
        enddo
        if(numdata < 2) return

        !Fill in the transform table
        allocate(ptt%table(2,numdata), stat = test)
        if(test /= 0)then
            write(*,'(/,a,/)') 'ERROR BUILDING TRANSTABLE, ALLOCATION FAILED'
            return
        endif
        j = 0
        do i = 1, nrow
            if(data_(colvr,i) >= tmin .and. data_(colvr,i) < tmax .and. &
               floor(data_(colct,i)+0.1) == ct)then
                j = j + 1
                ptt%table(1,j) = data_(colvr,i) + acorni(rng)*MACHP
                ptt%table(2,j) = 1.
            endif
        enddo

        !define size of the table
        ptt%n = numdata
        ptt%sumw = real(numdata,8)

        ! Compute nscores
        call nscore_transtable ( ptt, numq )

        build_transtable_by_category = .true.
        return
    end function build_transtable_by_category

    !Build a weighted transform table by category
    ! this routine is called from build_transtable so error checking is not needed
    logical function build_transtablew_by_category ( data_, colvr, colwt, colct, ct, numq, tmin, tmax, ptt )
        real*8, pointer :: data_(:,:) !a block of data for the table, must be ncol by nrows
        integer, intent (in) :: colvr, colwt, colct !variable and category columns
        integer, intent (in) :: ct !category of interest
        integer, intent (in) :: numq !number of quantiles
        real*8, intent (in) :: tmin, tmax !trimming limits
        type(transtable), intent(inout) :: ptt !pointer to the transform table to be built
        integer :: ncol, nrow
        integer :: i, j, test
        integer :: numdata
        real*8 :: sumweight

        build_transtablew_by_category = .false.
        ncol = ubound(data_,1)
        nrow = ubound(data_,2)
        numdata = 0
        sumweight = 0.

        !Count the data
        do i = 1, nrow
            if(data_(colvr,i) >= tmin .and. data_(colvr,i) < tmax .and. &
               floor(data_(colct,i)+0.1) == ct)then
                numdata = numdata + 1
            endif
        enddo
        if(numdata < 2) return

        !Fill in the transform table
        allocate(ptt%table(2,numdata), stat = test)
        if(test /= 0)then
            write(*,'(/,a,/)') 'ERROR BUILDING TRANSTABLE, ALLOCATION FAILED'
            return
        endif
        j = 0
        do i = 1, nrow
            if(data_(colvr,i) >= tmin .and. data_(colvr,i) < tmax .and. &
               floor(data_(colct,i)+0.1) == ct)then
                j = j + 1
                ptt%table(1,j) = data_(colvr,i) + acorni(rng)*MACHP
                ptt%table(2,j) = data_(colwt,i)
                if(ptt%table(2,j) < MACHP )then ! Check if the weight is numerically small
                    ptt%table(2,j) = MACHP
                endif
                sumweight = sumweight + ptt%table(2,j)
            endif
        enddo

        !define size of the table
        ptt%n = numdata
        ptt%sumw = sumweight

        ! Compute nscores
        call nscore_transtable ( ptt, numq )

        build_transtablew_by_category = .true.
        return
    end function build_transtablew_by_category

    ! Compute the normal scores for a transform table
    ! The table contains the number of data and total weight needed.
    ! If numq is passed in > 0, then the size of the table is
    ! cut back.
    subroutine nscore_transtable ( ptt, numq )
        type(transtable), target, intent(inout) :: ptt
        integer, intent (in) :: numq
        real*8 :: w, wf, cp0, cp1
        real*8, pointer :: z(:),q(:)
        integer :: i, j, n

        ! Sort data by value (this sort call sorts trtab by column 1)
        n = ptt%n
        i = qsort(ptt%table,1,n,(/1/))
        wf = 1. / ptt%sumw
        cp0 = 0.
        cp1 = 0.

        if(n > numq .and. numq > 0)then ! Resize for specified numq
            !Compute CDF
            do i = 1, n
                w =  wf * ptt%table(2,i)
                cp1 = cp1 + w
                ptt%table(2,i) = (cp1 + cp0) / 2.
                cp0 = cp1
            enddo

            !Generate a set of equally spaced normal values covering
            !the range of quantiles in the full size transform table
            !cp0 and cp1 are reused for miny and maxy and w as the diff
            allocate(rwork(2,numq))
            z => ptt%table(1,:)
            q => ptt%table(2,:)
            cp0 = snorm_inv(q(1))
            cp1 = snorm_inv(q(n))
            w = (cp1 - cp0) / real(numq-1,8)

            do i = 1, numq
                rwork(2,i) = snorm_cdf(cp0+real(i-1,8)*w)
                call locate(q,n,1,n,rwork(2,i),j)
                j = min(max(1,j),(n-1))
                rwork(1,i) = powint(q(j),q(j+1),z(j),z(j+1),rwork(2,i),1.D0)
                rwork(2,i) = snorm_inv ( rwork(2,i) )
            enddo
            rwork(1,numq) = ptt%table(1,n)

            deallocate(ptt%table)
            allocate(ptt%table(2,numq))
            ptt%n = numq
            ptt%table(1,:) = rwork(1,:)
            ptt%table(2,:) = rwork(2,:)

            !Add small number to max so there is room for breaking ties
            ptt%table(1,1) = ptt%mint
            ptt%table(1,numq) = ptt%table(1,numq) + MACHP

            nullify(z)
            nullify(q)
            deallocate(rwork)
        else
            ! Compute the cumulative probabilities and nscore values
            do i = 1, n
                w =  wf * ptt%table(2,i)
                cp1 = cp1 + w
                ptt%table(2,i) = snorm_inv ( (cp1 + cp0) / 2. )
                cp0 = cp1
            enddo

            !Add small number to max so there is room for breaking ties
            ptt%table(1,n) = ptt%table(1,n) + MACHP

        endif

        return
    end subroutine nscore_transtable

end module normalscore

max_stars_repo_name: haroldvelasquez/pygeostat

max_issues_repo_name: haroldvelasquez/pygeostat

max_forks_repo_name: haroldvelasquez/pygeostat