text:
      SUBROUTINE PLA2
C*****
C THIS ROUTINE IS THE SECOND FUNCTIONAL MODULE UNIQUE TO THE PIECE-WISE
C LINEAR ANALYSIS (PLA) RIGID FORMAT FOR THE DISPLACEMENT APPROACH.
C
C DMAP CALL...
C
C PLA2   DELTAUGV,DELTAPG,DELTAQG/UGV1,PGV1,QG1/V,N,PLACOUNT/ $
C
C CONCERNING DELTAUGV AND UGV1, THE ROUTINE WORKS AS FOLLOWS...
C DELTAUGV IS THE CURRENT INCREMENTAL DISPLACEMENT VECTOR IN THE PLA
C RIGID FORMAT LOOP AND UGV1 IS AN APPENDED FILE OF DISPLACEMENT VECTORS
C IF PLACOUNT .EQ. 1, THAT IS, THIS IS THE FIRST TIME PLA2 HAS BEEN
C CALLED IN THE PLA LOOP, THEN DELTAUGV IS COPIED ONTO UGV1.  IF
C PLACOUNT .GT. 1, THE PREVIOUS, OR LAST, DISPLACEMENT VECTOR IS READ
C INTO CORE FROM THE UGV1 DATA BLOCK, AND THE UGV1 FILE IS CLOSED WITH-
C OUT REWIND, THEN OPENED WITHOUT REWIND TO WRITE.  THE DELTAUGV VECTOR
C IS READ AN ELEMENT AT A TIME USING SUBROUTINE ZNTPKI AND ADDED TO
C THE VECTOR IN CORE.  THEN THE NEW DISPLACEMENT VECTOR IS WRITTEN ONTO
C THE UGV1 FILE.
C
C THEN THE PLA DMAP LOOP COUNTER PLACOUNT IS INCREMENTED.
C
C DELTAPG IS THE CURRENT INCREMENTAL LOAD VECTOR AND PGV1 IS THE
C CORRESPONDING MATRIX OF RUNNING SUM LOAD VECTORS.  PROCESSING IS
C SIMILAR TO THE ABOVE.  NOTE THAT NEITHER DATA BLOCK, LIKE THE TWO
C DISCUSSED ABOVE, CAN BE PURGED.
C
C DELTAQG IS THE CURRENT INCREMENTAL VECTOR OF SINGLE POINT CONSTRAINT
C FORCES AND QG1 IS THE APPENDED FILE OF VECTORS OF SPCF.  THESE TWO
C DATA BLOCKS ARE PROCESSED IDENTICALLY TO DELTAUGV AND UGV1 EXCECT
C THAT NO FATAL ERROR EXISTS IF ONE OR THE OTHER HAS BEEN PURGED.
C*****
C
      INTEGER
     1                   BUFSZ
     2,                  BUFFR1             ,BUFFR2
     3,                  OFILE              ,PLACNT
     4,                  EOR                ,CLSRW
     5,                  CLSNRW             ,OUTRW
     6,                  EOL                ,TYPEA
     7,                  TYPEB              ,OUTNRW
      INTEGER INBLK(15),OUBLK(15)
C
      DIMENSION
     1                   NAME(2)            ,DUMMY(2)
     2,                  MCB(7)
      COMMON /BLANK/PLACNT
      COMMON   /SYSTEM/  BUFSZ
      COMMON   /ZZZZZZ /  Z(1)
      COMMON   /ZNTPKX/
     1                   A(4)               ,INDEX
     2,                  EOL                ,IDUMMY
      COMMON   /PACKX /
     1                   TYPEA              ,TYPEB
     2,                  IPACK              ,JPACK
     3,                  INCPK
      COMMON   /UNPAKX/
     1                   IUNPKB             ,IUNPK
     2,                  JUNPK              ,INCUPK
C
      DATA               NAME /4HPLA2,4H    /
      DATA               INRW,OUTRW,OUTNRW,CLSRW,CLSNRW,EOR/0,1,3,1,2,1/
C
C INITIALIZE
C
      IZMAX = KORSZ (Z)
      BUFFR1 = IZMAX - BUFSZ
      BUFFR2 = BUFFR1 - BUFSZ
      LEFT = BUFFR2 - 1
      ILOOP = 1
      IFILE = 101
      OFILE = 201
C
C OPEN INPUT FILE TO READ AND OUTPUT FILE TO WRITE (IF PLACNT .EQ. 1)
C OR TO READ (IF PLACNT .GT. 1)
C
   10 JFILE = IFILE
      INBLK(1) = IFILE
      OUBLK(1) = OFILE
      DO 15 I = 2,7
   15 MCB(I) = 0
      MCB(1) = OFILE
      IF (PLACNT .EQ. 1) MCB(1) = IFILE
      CALL RDTRL (MCB)
      CALL OPEN(*100,IFILE,Z(BUFFR1),INRW)
      CALL FWDREC(*9020,IFILE)
      IOPT = INRW
      IF (PLACNT .EQ. 1) IOPT = OUTRW
      CALL OPEN(*110,OFILE,Z(BUFFR2),IOPT)
      IF (PLACNT .NE. 1) GO TO 30
C
C THIS IS THE FIRST TIME THROUGH THE PLA LOOP.  COPY THE VECTOR ON THE
C INPUT FILE ONTO THE OUTPUT FILE.
C
      CALL FNAME (OFILE,DUMMY)
      CALL WRITE (OFILE,DUMMY,2,EOR)
      CALL CPYSTR(INBLK,OUBLK,0,0)
      CALL CLOSE (OFILE,CLSRW)
      CALL CLOSE(IFILE,CLSRW)
      GO TO 70
C
C THIS IS NOT THE FIRST PASS
C
   30 JFILE = OFILE
      CALL FWDREC(*9020,OFILE)
      NRECS = PLACNT - 2
      IF (NRECS .LE. 0) GO TO 50
      DO 40 I = 1,NRECS
   40 CALL FWDREC(*9020,OFILE)
   50 MCB(1) = OFILE
      CALL RDTRL (MCB)
      MCB(6) = 0
      MCB(7) = 0
      IF (LEFT .LT. MCB(3)) CALL MESAGE (-8,0,NAME)
      IUNPKB = 1
      IUNPK  = 1
      JUNPK  = MCB(3)
      INCUPK = 1
      CALL UNPACK(*9030,OFILE,Z)
      CALL CLOSE (OFILE,CLSNRW)
      CALL OPEN(*9010,OFILE,Z(BUFFR2),OUTNRW)
C
C READ THE INCREMENTAL VECTOR.  INTPK INITIALIZES AND ZNTPKI RETURNS
C ONE ELEMENT AT A TIME
C
      KTYPE = 1
      CALL INTPK(*9040,IFILE,0,KTYPE,0)
C
C READ AND ADD LOOP.
C
   60 CALL ZNTPKI
      Z(INDEX) = Z(INDEX) + A(1)
      IF (EOL .EQ. 0) GO TO 60
C
C ADDITION HAS BEEN COMPLETED
C
      CALL CLOSE (IFILE,CLSRW)
C
C WRITE VECTOR ON OUTPUT FILE IN PACKED FORMAT.
C
      TYPEA = 1
      TYPEB = 1
      IPACK = 1
      JPACK = MCB(3)
      INCPK = 1
      CALL PACK(Z,OFILE,MCB)
      CALL CLOSE (OFILE,CLSRW)
C
C WRITE TRAILER
C
   70 MCB(1) = OFILE
      CALL WRTTRL (MCB)
   90 ILOOP = ILOOP + 1
      IF (ILOOP .GT. 3) GO TO 200
      IFILE = IFILE + 1
      OFILE = OFILE + 1
      GO TO 10
  100 IF (ILOOP .EQ. 1  .OR.  ILOOP .EQ. 2) CALL MESAGE (-30,127,IFILE)
      GO TO 90
  110 IF (ILOOP .EQ. 1  .OR.  ILOOP .EQ. 2) CALL MESAGE (-30,128,OFILE)
      GO TO 90
C
C INCREMENT THE PLA DMAP LOOP COUNTER
C
  200 PLACNT = PLACNT + 1
      RETURN
C
C FATAL ERRORS
C
 9010 CALL MESAGE (-1,JFILE,NAME)
 9020 CALL MESAGE (-2,JFILE,NAME)
 9030 CALL MESAGE (-30,129,ILOOP+200)
 9040 CALL MESAGE (-30,130,ILOOP+100)
      RETURN
      END


text:
module problem_io_module

  implicit none

  ! For determining if we are the I/O processor.

  logical, save :: ioproc

  ! Probin file

  character (len=:), allocatable, save :: probin

contains

  subroutine initialize_io(name, namlen)

    implicit none

    integer :: namlen, i
    integer :: name(namlen)

    integer :: is_ioprocessor

    ! Build "probin" filename -- the name of the file containing the fortin namelist.

    allocate(character(len=namlen) :: probin)
    do i = 1, namlen
       probin(i:i) = char(name(i))
    enddo

    ! Determine whether we are the I/O procoessor.

    call bl_pd_is_ioproc(is_ioprocessor)

    ioproc = is_ioprocessor .ne. 0

  end subroutine initialize_io

end module problem_io_module

max_stars_repo_name: yingtchen/Castro

max_issues_repo_name: yingtchen/Castro

max_forks_repo_name: yingtchen/Castro


text:
program problem78
  use fmzm
  implicit none
  integer            :: i
  type(im)           :: temp,divisor,neg_one,zero,one,two
  integer, parameter :: memo_size=100000
  type(im)           :: memo(0:memo_size)

  neg_one=to_im(-1)
  zero=to_im(0)
  one=to_im(1)
  two=to_im(2)
  divisor=to_im(1000000)

  ! https://oeis.org/A000041
  do i=1,memo_size
     memo(i)=to_im(-1)
  end do
  memo(0)=one

  do i=1,memo_size
     temp=partitions(i)
     if (mod(temp,divisor)==zero) then
        print *, i
        exit
     end if
  end do

contains

  pure function plus_minus(n)
    implicit none
    integer, intent(in) :: n
    integer             :: plus_minus

    plus_minus=(-1)**((n-1)/2)

  end function plus_minus

  ! https://oeis.org/A001318
  function pentagonal(n)
    implicit none
    integer, intent(in) :: n
    integer             :: m,pentagonal

    if (mod(n,2)==1) then
       m=1
    else
       m=-1
    end if
    m=m*(n+1)/2
    pentagonal=m*(3*m-1)/2

  end function pentagonal

  ! https://en.wikipedia.org/wiki/Partition_(number_theory)#Generating_function
  function partitions(n)
    use fmzm
    implicit none
    integer, intent(in) :: n
    integer             :: i
    type(im)            :: partitions

    if (n<0) then
       partitions=zero
       return
    else if (n==0) then
       partitions=one
       return
    else if (memo(n)>neg_one) then
       partitions=memo(n)
       return
    else
       partitions=zero
       i=1
       do while (pentagonal(i)<=n)
          partitions=partitions+to_im(plus_minus(i))*memo(n-pentagonal(i))
          i=i+1
       end do
       memo(n)=partitions
       return
    end if

  end function partitions

end program problem78

max_stars_repo_name: jamesmcclain/Dioskouroi

max_issues_repo_name: jamesmcclain/Dioskouroi

max_forks_repo_name: jamesmcclain/Dioskouroi


text:
Tapas are essentially appetizers in Spanish cuisine, including both cold and hot foods.  They often include cheese or fruit platters, fried croquettes, grilled meats, and so forth.  Theyre well suited for groups of any size; a meal typically consists of 1.52 tapas per person (though this may vary by location), which are typically shared.

Sadly, Davis does not have any dedicated tapas restaurants.   Options in the area include:

  Ficelle in Winters
  Tapa the World in Midtown Sacramento
  Aioli Bodega Espanola in Sacramento (the Davis location closed)

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
    !Proposal density
    !Using the covariance matrix to give the proposal directions typically
    !significantly increases the acceptance rate and gives faster movement
    !around parameter space.
    !We generate a random basis in the eigenvectors, then cycle through them
    !proposing changes to each, then generate a new random basis
    !The distance proposal in the random direction is given by a two-D Gaussian
    !radial function mixed with an exponential, which is quite robust to wrong width estimates
    !See http://cosmologist.info/notes/cosmomc.pdf

    module propose
    use RandUtils
    use settings
    use GeneralTypes
    implicit none

    Type IndexCycler
        integer n
        integer :: loopix=0
    end Type IndexCycler

    Type, extends(IndexCycler) :: CyclicIndexRandomizer
        integer, allocatable :: indices(:)
    contains
    procedure :: Next
    end Type CyclicIndexRandomizer

    Type, extends(IndexCycler) :: RandDirectionProposer
        integer :: propose_count = 0
        real(mcp), allocatable, dimension(:,:) :: R
    contains
    procedure :: ProposeVec
    procedure :: Propose_r
    end type RandDirectionProposer

    Type, extends(RandDirectionProposer) :: BlockProposer
        integer :: block_start
        integer, allocatable :: used_param_indices(:)
        integer, allocatable :: params_changed(:), used_params_changed(:)
        real(mcp), allocatable :: mapping_matrix(:,:)
    contains
    procedure :: UpdateParams
    end Type BlockProposer

    Type, extends(TCheckpointable) :: TCovmatProposer
        real(mcp), allocatable :: propose_matrix(:,:)
    contains
    procedure :: SetCovariance => TCovmatProposer_SetCovariance
    procedure :: SaveState => TCovmatProposer_SaveState
    procedure :: LoadState => TCovmatProposer_LoadState
    end type

    Type, extends(TCovmatProposer) :: BlockedProposer
        integer :: nBlocks
        integer, allocatable :: indices(:), ProposerForIndex(:)
        Type(CyclicIndexRandomizer) :: Slow, Fast, All
        Type(BlockProposer), allocatable :: Proposer(:)
        integer :: oversample_fast = 1
        real(mcp) :: propose_scale = 2.4_mcp
        integer, private  :: fast_ix = 0
    contains
    procedure :: Init
    procedure :: SetCovariance => BlockedProposer_SetCovariance
    procedure :: GetBlockProposal
    procedure :: GetProposal
    procedure :: GetProposalSlow
    procedure :: GetProposalFast
    procedure :: GetProposalFastDelta
    end Type BlockedProposer

    logical :: propose_rand_directions = .true.

    contains

    function Next(this)
    Class(CyclicIndexRandomizer) :: this
    integer Next

    this%loopix = mod(this%loopix, this%n) + 1
    if (this%loopix == 1) then
        if (.not. allocated(this%indices)) allocate(this%indices(this%n))
        call randIndices(this%indices, this%n, this%n)
    end if
    Next = this%indices(this%loopix)

    end function Next

    subroutine RotMatrix(M,n)
    integer, intent(in) :: n
    real(mcp) M(n,n)
    integer i

    if (propose_rand_directions .and. n > 1) then
        call RandRotation(M, n)
    else
        M = 0
        do i = 1, n
            M(i,i) = sign(1._mcp, real(ranmar(),mcp)-0.5_mcp)
        end do
    end if

    end subroutine RotMatrix

    function ProposeVec(this, wid) result(vec)
    Class(RandDirectionProposer) :: this
    real(mcp) :: vec(this%n)
    real(mcp) :: wid

    if (mod(this%loopix,this%n)==0) then
        !Get a new random rotation
        if (.not. allocated(this%R)) allocate(this%R(this%n,this%n))
        call RotMatrix(this%R, this%n)
        this%loopix = 0
    end if
    this%loopix = this%loopix + 1
    this%propose_count= this%propose_count + 1
    vec = this%R(:,this%loopix) * ( this%Propose_r() * wid )

    end function ProposeVec

    function Propose_r(this) result(r_fac)
    !distance proposal function (scale is 1)
    Class(RandDirectionProposer), intent(in) :: this
    integer i,n
    real(mcp) r_fac

    if (ranmar() < 0.33d0) then
        r_fac = randexp1()
    else
        n = min(this%n, 2)
        r_fac = 0
        do i = 1, n
            r_fac = r_fac + Gaussian1()**2
        end do
        r_fac = sqrt( r_fac / n )
    end if

    end function Propose_r

    subroutine UpdateParams(this, P, vec)
    Class(BlockProposer) :: this
    real(mcp) :: P(:)
    real(mcp):: vec(this%n)

    P(this%params_changed) =  P(this%params_changed) + matmul(this%mapping_matrix, vec)

    end subroutine UpdateParams

    subroutine Init(this, parameter_blocks, slow_block_max, oversample_fast, propose_scale)
    !slow_block_max determines which parameter blocks are grouped together as being 'slow'
    Class(BlockedProposer), target :: this
    type(int_arr) :: parameter_blocks(:)
    integer, intent(in) :: slow_block_max
    integer used_blocks(size(parameter_blocks))
    integer, intent(in), optional :: oversample_fast
    real(mcp), intent(in), optional :: propose_scale
    integer i, ix, n, np
    Type(BlockProposer), pointer :: BP

    if (present(oversample_fast)) then
        this%oversample_fast = oversample_fast
    endif
    if (present(propose_scale)) then
        this%propose_scale = propose_scale
    end if
    this%nBlocks = size(parameter_blocks)
    n=0
    this%All%n=0
    this%Slow%n=0
    do i=1, this%nBlocks
        np = size(parameter_blocks(i)%P)
        if (np >0) then
            this%All%n = this%All%n+ np
            if (i <= slow_block_max) this%Slow%n= this%Slow%n + np
            n=n+1
            used_blocks(n)=i
        end if
    end do
    this%Fast%n = this%All%n - this%Slow%n
    this%nBlocks = n
    allocate(this%Proposer(this%nBlocks))
    allocate(this%indices(this%All%n))
    allocate(this%ProposerForIndex(this%All%n))
    this%indices=0
    ix = 1
    do i=1, this%nBlocks
        !loop from slow to fast
        BP => this%Proposer(i)
        BP%block_start = ix
        BP%n = size(parameter_blocks(used_blocks(i))%P)
        allocate(BP%used_param_indices(BP%n))
        BP%used_param_indices = parameter_blocks(used_blocks(i))%P
        this%indices(ix:ix+BP%n-1) = BP%used_param_indices
        this%ProposerForIndex(ix:ix+BP%n-1) = i
        ix = ix + BP%n
    end do
    if (any(this%indices==0)) stop 'DecomposeCovariance: not all used parameters blocked'
    do i=1, this%nBlocks
        BP => this%Proposer(i)
        allocate(BP%used_params_changed(this%All%n-BP%block_start+1))
        allocate(BP%params_changed(this%All%n-BP%block_start+1))
        BP%used_params_changed = this%indices(BP%block_start:this%All%n)
        BP%params_changed = params_used(BP%used_params_changed)
    end do

    end subroutine Init

    subroutine BlockedProposer_SetCovariance(this, propose_matrix)
    !take covariance of used parameters (propose_matrix), and construct orthonormal parmeters
    !where orthonormal parameters are grouped in blocks by speed, so changes in slowest block
    !changes slow and fast parameters, but changes in the fastest block only changes fast parameters
    use MatrixUtils
    Class(BlockedProposer), target :: this
    real(mcp), intent(in) :: propose_matrix(:,:)
    real(mcp) corr(num_params_used,num_params_used)
    real(mcp) :: sigmas(num_params_used), L(num_params_used,num_params_used)
    integer i, j
    Type(BlockProposer), pointer :: BP

    if (.not. allocated(this%propose_matrix)) then
        allocate(this%propose_matrix, source=propose_matrix)
    else
        this%propose_matrix=propose_matrix
    end if
    do i = 1, num_params_used
        sigmas(i) = sqrt(propose_matrix(i,i))
        corr(i,:) = propose_matrix(i,:) / sigmas(i)
    end do
    do i = 1, num_params_used
        corr(:,i) = corr(:,i) / sigmas(i)
    end do
    L = corr(this%indices,this%indices)
    call Matrix_Cholesky(L,zeroed=.true.)
    do i=1, this%nBlocks
        BP => this%Proposer(i)
        if (.not. allocated(BP%mapping_matrix)) allocate(BP%mapping_matrix(size(BP%used_params_changed), BP%n))
        do j = 1, size(BP%used_params_changed)
            BP%mapping_matrix(j,:)  =  sigmas(BP%used_params_changed(j)) * &
            L(BP%block_start+j-1,BP%block_start:BP%block_start+BP%n-1)
        end do
    end do
    end subroutine BlockedProposer_SetCovariance

    subroutine GetBlockProposal(this, P, i)
    class(BlockedProposer) :: this
    integer, intent(in) :: i !i must be first in gfortran
    real(mcp) :: P(:)

    call this%Proposer(i)%UpdateParams(P, this%Proposer(i)%Proposevec(this%propose_scale))

    end subroutine GetBlockProposal

    subroutine GetProposal(this, P)
    class(BlockedProposer) :: this
    real(mcp) :: P(:)

    if (this%fast_ix/=0) then
        call  GetProposalFast(this, P)
        this%fast_ix=this%fast_ix-1
    else
        if (this%All%Next() > this%Slow%n) then
            call  GetProposalFast(this, P)
            this%fast_ix = this%oversample_fast-1
        else
            call  GetProposalSlow(this, P)
        end if
    end if

    end subroutine GetProposal

    subroutine GetProposalSlow(this, P)
    class(BlockedProposer) :: this
    real(mcp) :: P(:)

    call this%GetBlockProposal(P, this%ProposerForIndex(this%Slow%Next()))

    end subroutine GetProposalSlow

    subroutine GetProposalFast(this, P)
    class(BlockedProposer) :: this
    real(mcp) :: P(:)

    call this%GetBlockProposal(P, this%ProposerForIndex(this%Slow%n + this%Fast%Next()))

    end subroutine GetProposalFast

    subroutine GetProposalFastDelta(this, P)
    class(BlockedProposer) :: this
    real(mcp) :: P(:)

    P=0
    call this%GetProposalFast(P)

    end subroutine GetProposalFastDelta

    subroutine TCovmatProposer_SetCovariance(this, propose_matrix)
    class(TCovmatProposer), target :: this
    real(mcp), intent(in) :: propose_matrix(:,:)
    !only inherited otherwise do nothing
    end subroutine TCovmatProposer_SetCovariance

    subroutine TCovmatProposer_LoadState(this,F)
    class(TCovmatProposer) :: this
    class(TFileStream) :: F
    real(mcp) :: cov(num_params_used,num_params_used)

    if (.not. F%ReadItem(cov)) call MpiStop('error reading proposal covariance')
    call this%SetCovariance(cov)

    end subroutine TCovmatProposer_LoadState

    subroutine TCovmatProposer_SaveState(this,F)
    class(TCovmatProposer) :: this
    class(TFileStream) :: F

    call F%Write(this%propose_matrix)

    end subroutine TCovmatProposer_SaveState

    end module propose

max_stars_repo_name: sjoudaki/CosmoJBD

max_issues_repo_name: sjoudaki/CosmoJBD

max_forks_repo_name: sjoudaki/CosmoJBD


text:
MODULE alarrays

  use globvars

!***********************************************************************************!
!*
!*         Allocation/Deallocation Module
!*
!*   Contains subroutines for:
!*
!*      1) Allocating & initialising the basis set
!*      2) Allocating & initialising a single basis function
!*      3) Deallocating a single basis function
!*      4) Deallocating the basis set
!*      5) Allocating and initialising the multi-configurational nbf x nbf
!*           Hamiltonian matrix
!*      6) Allocating and initialising the single configurational npes x npes
!*           Hamiltonian matrices
!*      7) Deallocating the multi-configurational Hamiltonian matrix
!*      8) Deallocating the single-configurational Hamiltonian matrices
!*
!***********************************************************************************!

contains

!***********************************************************************************!
!               Array Allocation/Deallocation Subroutines
!***********************************************************************************!

  subroutine allocbs(bs, nbf)

    ! Allocates the top level of the variable of the basisfn defined type

    implicit none

    type(basisfn),dimension(:),allocatable, intent(inout)::bs
    integer, intent (in) :: nbf

    integer:: j, ierr

    if (errorflag .ne. 0) return

    ierr=0
    if (allocated(bs).eqv..false.) then
      allocate (bs(nbf), stat=ierr) !allocate the top level of the defined type array
      if (ierr/=0) then
        write(0,"(a,i0)") "Error in basis set allocation. ierr had value ", ierr
        errorflag=1
        return
      end if
    end if

    do j=1,size(bs)
      call allocbf(bs(j))! this allocates each bf within bs
    end do

    return

  end subroutine allocbs

!------------------------------------------------------------------------------------

  subroutine allocbf(bf)

    ! Allocates the elements of the variable of the basisfn defined type

    implicit none

    type(basisfn),intent(inout)::bf

    integer :: ierr

    if (errorflag .ne. 0) return

    ierr=0

    ! Each sub-array must be allocated separately.
    ! Allocation continues only if no errors occur in previous allocations

    allocate (bf%z(ndim), stat=ierr)
    if (ierr==0) allocate (bf%d_pes(npes), stat=ierr)
    if (ierr==0) allocate (bf%s_pes(npes), stat=ierr)
    if (ierr==0) allocate (bf%a_pes(npes), stat=ierr)
    if (ierr/=0) then
      write(0,"(a,i0)") "Error in basis function allocation. ierr had value ", ierr
      errorflag=1
      return
    end if

    bf%z    (1:ndim) = (0.0d0,0.0d0)
    bf%d_pes(1:npes) = (0.0d0,0.0d0)
    bf%a_pes(1:npes) = (0.0d0,0.0d0)
    bf%s_pes(1:npes) =  0.0d0
    bf%D_big         = (0.0d0,0.0d0)

    return

  end subroutine allocbf

!------------------------------------------------------------------------------------

  subroutine deallocbf(bf)

    ! Deallocates the elements of the variable of the basisfn defined type

    implicit none

    type(basisfn),intent(inout)::bf

    integer :: ierr

    if (errorflag .ne. 0) return

    ierr=0

    ! Each sub-array must be deallocated separately.
    ! Deallocation continues only if no errors occur in previous allocations

    deallocate (bf%z, stat=ierr)
    if (ierr==0) deallocate (bf%d_pes, stat=ierr)
    if (ierr==0) deallocate (bf%s_pes, stat=ierr)
    if (ierr==0) deallocate (bf%a_pes, stat=ierr)
    if (ierr/=0) then
      write(0,"(a,i0)") "Error in basis function deallocation. ierr had value ", ierr
      errorflag=1
      return
    end if

    return

  end subroutine deallocbf

!------------------------------------------------------------------------------------

  subroutine deallocbs(bs)

    ! Deallocates the top level of the variable of the basisfn defined type

    implicit none

    type(basisfn),dimension(:), allocatable, intent(inout)::bs

    integer:: j, ierr

    if (errorflag .ne. 0) return

    ierr=0

    do j=1,size(bs)
      call deallocbf(bs(j))! this deallocates each bf within bs
    end do

    ! Deallocation is carried out bottom up.
    ! Top level deallocation only occurs if lower level deallocation was successful

    if (errorflag==0) deallocate (bs, stat=ierr)
    if (ierr/=0) then
      write(0,"(a,i0)") "Error in basis set deallocation. ierr had value ", ierr
      errorflag=1
      return
    end if

    return

  end subroutine deallocbs

!------------------------------------------------------------------------------------

  subroutine allocham(H, nbf)

    ! Allocates the top level of the variable of the hamiltonian defined type

    implicit none

    type(hamiltonian),dimension(:,:),allocatable,intent (inout)::H
    integer, intent (in) :: nbf

    integer::k, j, ierr

    if (errorflag .ne. 0) return

    ierr=0

    ! Allocate the top level of the array before calling the subroutine to allocate
    ! the different sub-arrays

    allocate (H(nbf, nbf), stat=ierr)
    do k=1,nbf
      do j=1,nbf
        call allochamjk(H(j,k))
      end do
    end do
    if (ierr/=0) then
      write(0,"(2a,i0)") "Error in hamiltonian matrix allocation.", &
                          " ierr had value ", ierr
      errorflag=1
      return
    end if

  end subroutine allocham

!------------------------------------------------------------------------------------

  subroutine allochamjk(Hsing)

    ! Allocates the elements of the variable of the hamiltonian defined type

    implicit none

    type(hamiltonian),intent(inout) :: Hsing

    integer :: ierr

    if (errorflag .ne. 0) return

    ierr = 0

    ! Each sub-array has to be allocated separately

    allocate (Hsing%Hjk(npes, npes), stat=ierr)
    if (ierr/=0) then
      write(0,"(a,i0)") "Error in hamiltonian element allocation.",&
                          " ierr had value ", ierr
      errorflag=1
      return
    end if
    Hsing%Hjk = (0.0d0,0.0d0)

    return

  end subroutine allochamjk

!------------------------------------------------------------------------------------

  subroutine deallocham(H)

    ! Deallocates the top level of the variable of the hamiltonian defined type

    implicit none

    type(hamiltonian),dimension(:,:),allocatable,intent (inout)::H

    integer::k, j, ierr

    if (errorflag .ne. 0) return

    ierr=0

    ! Deallocate the top level of the array before calling the subroutine to
    ! deallocate the different sub-arrays

    do k=1,size(H,2)
      do j=1,size(H,1)
        call deallochamjk(H(j,k))
      end do
    end do
    deallocate (H, stat=ierr)
    if (ierr/=0) then
      write(0,"(a,i0)") "Error in hamiltonian matrix deallocation.",&
                          " ierr had value ", ierr
      errorflag=1
      return
    end if

  end subroutine deallocham

!------------------------------------------------------------------------------------

  subroutine deallochamjk(Hsing)

    ! Deallocates the elements of the variable of the hamiltonian defined type

    implicit none

    type(hamiltonian),intent(inout) :: Hsing

    integer :: ierr

    if (errorflag .ne. 0) return

    ierr = 0

    ! Each sub-array is deallocated separately

    deallocate (Hsing%Hjk, stat=ierr)
    if (ierr/=0) then
      write(0,"(a,i0)") "Error in hamiltonian element deallocation.",&
                          " ierr had value ", ierr
      errorflag=1
      return
    end if

    return

  end subroutine deallochamjk

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

END MODULE alarrays

max_stars_repo_name: Leeds-MCE/MCE

max_issues_repo_name: Leeds-MCE/MCE

max_forks_repo_name: Leeds-MCE/MCE


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program 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 General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

!
!  Subroutine DELAY (NSEC)
!  This program pauses execution for NSEC seconds, unless the user hits
!  any key, in which case execution continues immediately.
!
!  Program written by Hans Los
!  Version 1.0
!  Released july 1989.
!  Part of the INTERACT.LIB.
!
!
      SUBROUTINE DELAY (NSEC)
      INTEGER*2 HOUR, MINUTE, SEC, SEC100
      INTEGER*4 SECMID, SECEND, NSEC, KEY
!
      CALL GETTIM (HOUR, MINUTE, SEC, SEC100)
      SECMID = HOUR * 3600 + MINUTE * 60 + SEC
      SECEND = SECMID + NSEC
10    CONTINUE
      CALL CHKKEY (KEY)
      IF (KEY .NE. 0) GO TO 20
      CALL GETTIM (HOUR, MINUTE, SEC, SEC100)
      SECMID = HOUR * 3600 + MINUTE * 60 + SEC
      IF (SECMID .LT. SECEND) GO TO 10
20    CONTINUE
      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
A saloon was a type of bars bar popular in the American West prior to the Prohibition. There were many saloons in Davisville up until 1911 when they were finally banned due largely to the efforts of the Womens Christian Temperance Union. Many of the saloons had bad reputations and were associated with violent crime in particular. Two saloonkeepers were shot in the head by angry customers on separate occasions, but both somehow survived and made full recoveries.

  Chinese Saloon  near the 300 block of G Street Olive Street
  Depot Saloon  owned by Henry Anderson
  Fosters Saloon  owned by T. Foster
  Haly & Mayos Saloon  on 1st Street Front Street, near B Street Ash Street
  McLellan  on G Street Olive Street
  Ranchers Saloon  on G Street Olive Street, owned by John Caulfield
  Shillabers Saloon  on G Street Olive Street, owned by Oliver E. Dickerson
  Schultzs Saloon  near railroad, owned by Frederick Schultz
  http://www.woodlandwiki.org/The_Burger_Saloon The Burger Saloon  on Main street in Woodland, CA
  Webers Yolo Saloon  first true saloon in Davis history, on corner of 1st Street Front Street and G Street Olive Street, owned by George Weber
  Yolo House  a hotel but also the first place in Davis history to serve alcohol, on 1st Street Front Street and A Street Beech Street, operated by William Fredrich Dresbach William Dresbach

For more information on the history of drinking in Davisville, visit the town history entry. Also, see alcohol for a more modern take on catching a buzz.

There was once a Hair Saloon in Davis.  Unfortunately, it didnt serve alcohol.


text:
A dazzling array of Booya!

I have a ferocious apatite for all that is rad and a tendency to misspell everything.

Went to UC Davis UCD got my degree now I teach special ed.  I recently got into the JET program and I am leaving to Japan August 1st.

I have been a personal trainer at the ARC, an Outdoor guide and a swim instructor.  As much as I detest the term,
I am a foodie, and I love trying new food where I go, so most of my comments will/ have been about food.


text:
      subroutine gencol(x,xjac,xmin,emit,r)
      implicit none
c---Generate an x value and store it for later retrieval
      integer emit,lemit
      double precision x,xjac,xmin,xl,xljac,xlmin,r
      save xl,xljac,xlmin,lemit
      x=1d0-(1d0-xmin)*abs(1d0-2d0*r)
      xjac=2d0*(1d0-xmin)
      xl=x
      xljac=xjac
      xlmin=xmin
      lemit=emit
      return

      entry getcol(x,xjac,xmin,emit)
c---return the same values as last time
      x=xl
      xjac=xljac
      xmin=xlmin
      emit=lemit
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
       SUBROUTINE GATHED(SCR,LSCR,ISCR)
C                              PROCESS GATHER
C                              ------- ------
C
C  DOCUMENT DATE: 15 March 1994
C
C    A GATHER IS A COLLECTION OR REARRANGEMENT OF TRACES ACCORDING TO SOME
C  CRITERIA.  PROCESS GATHER COLLECTS OR REARRANGES THE INPUT TRACES ACCORDING
C  TO THE REFLECTION POINT (RP) NUMBER CALCULATED BY PROCESS GEOM.  THE GEOM
C  PARAMETERS MAY BE MANIPULATED BY THE USER TO GATHER THE INPUT TRACES ACCORDING
C  TO ANY CRITERIA BY FUDGING THE GEOM PARAMETERS.  A CONSTANT OFFSET GATHER
C  OF A UNIFORM MARINE LINE MAY BE MADE BY OMITTING TRACES VIA PROCESS INPUT.
C      PROCESS GATHER SORTS EACH GATHER BY THE ABSOLUTE VALUE OF THE SHOT-
C  RECEIVER DISTANCE, SO THAT THE SHORTEST RANGE TRACE IS FIRST WITHIN THE
C  GATHER.  EACH GATHER IS IS TERMINATED BY SETTING A SPECIAL FLAG IN THE
C  TRACE HEADER.  A GATHER RECORD IS THE COLLECTION OF ALL THESE TRACES.
C     SEE PROCESS GEOM FOR THE METHOD OF CALCULATING RP NUMBERS.
c     Gather creates a disk file to store the partial gathers while
c  the data is being read.  Gather assumes that the geometry of the data
c  does not skip around very much. i.e. the geometry doesn't go backwards
c  nor does it skip more than a cable length forward.  The temporary
c  disk file can hold MAXRPS rps (preset to 5 plus the number of traces
c  per shot from the SEGY binary header), with each each rp able to hold
c  a maximum of MAXTRS traces (also preset to the number of traces per
c  shot in the SEGY binary header), with each traces having a maximum
c  of NWRDS samples, with each sample being 4 bytes long (except on the
c  Cray).  The temporary disk file size will be:
c      maxtrs * maxrps * (nwrds+240) * 4
c  The preset values of maxtrs, maxrps, and nwrds are designed for the
c  marine geometry of advancing .5 groups between shots (96 cdp from
c  a 96 trace streamer).
C      A NULL SET OF GATHER PARAMETERS MUST BE GIVEN EVEN IF ALL THE
C  PARAMETERS ARE THE PRESETS.  E.G.  GATHER
C                                          END
C                                     END
C
C  THE PARAMETER DICTIONARY
C  --- --------- ----------
C  FRP    - THE FIRST RP NUMBER TO GATHER.  TRACES WITH RP NUMBERS LESS THAN FRP
C           ARE NOT GATHERED.  RP NUMBERS ARE CALCULATED BY PROCESS GEOM.
C           THIS PARAMETER SHOULD BE USED WHEN RESTARTING A PARTIALLY
C           COMPLETED GATHER RUN.  THE USE OF FRP WILL CAUSE PROCESS GATHER TO
C           REWRITE THE PROCESS INPUT PARAMETERS SO THAT THE UNNEEDED INPUT
C           WILL NOT BE PROCESSED.  BOTH THE PROCESS INPUT AND PROCESS GEOM
C           PARAMETERS MUST BE GIVEN PRIOR TO THE PROCESS GATHER PARAMETERS.
c           FRP only works with PROCESS INPUT.
C           PRESET=RP NUMBER OF THE FIRST TRACE.
C  RPINC  - THE INCREMENT OF RP NUMBERS BETWEEN THE RPS TO GATHER.  THE ONLY
C           TRACES GATHERED WILL HAVE BIN NUMBERS FRP, FRP+RPINC, FRP+2*RPINC,
C           FRP+3*RPINC, . . . . ETC.
C           PRESET=1.
C  NWRDS  - THE LARGEST NUMBER OF SAMPLES PER TRACE IN THE JOB.  THIS SHOULD BE
C           TRACE LENGTH PLUS THE TRACE HEADER LENGTH.
C           PRESET=FROM FIRST INPUT TRACE.
C  MAXRPS - THE MAXIMUM NUMBER OF BINS (OR RP'S) THAT ARE NEEDED ON
C           THE DISC AT ANY ONE TIME.  IN MARINE WORK THE NUMBER OF TRACES
C           PER SHOT WILL SUFFICE SINCE NO TWO UNGATHERED TRACES WITH THE SAME
C           RP NUMBER ARE MORE THAN A CABLE LENGTH AWAY.
C           PRESET=THE NUMBER OF TRACES FROM THE SEGY BINARY HEADER PLUS 5
C  MAXTRS - THE MAXIMUM NUMBER OF TRACES ANY ONE GATHER CAN HAVE.  IN RP GATHERS
C           THIS IS THE MAXIMUM CDP ALLOWED.
C           PRESET=THE NUMBER OF TRACES FROM THE SEGY BINARY HEADER.
C  MINTRS - THE MINIMUM NUMBER OF TRACES EACH GATHER CAN HAVE.  IF MINTRS=0 AND
C           NO INPUT TRACES CONTRIBUTE TO A GIVEN GATHER, THAT GATHER WILL NOT
C           BE OUTPUT.
C           PRESET=1    e.g. mintrs 24
C
C  WRITTEN AND COPYRIGHTED BY:
C   PAUL HENKART, SCRIPPS INSTITUTION OF OCEANOGRAPHY, MARCH 1980
C
C
C   MODIFIED ON 3 MAY 1981 FOR THE RESTART CAPABILITY.  P.C.H.
c   modified on 24 May 1989 so the restart will work under Version 2.0
c     ( open the new input parameter file correctly and read/write
c      the correct items for the Version 2 process input parameters).
c   mod 3 Dec 1993.  Preset nwrds to 0 and let the execute determine the
c     length from the first GOOD trace it gets.
c   mod 20 Apr 96 - Skip the restart frp if geom isn't there and the cdp
c      number is is in the trace already (like geom was done in a
c      different run)
c  mod 18 Nov 98 - Allow frp without process input, in case the first
c      trace read is not the first trace of the subsurface.
c  mod March 1999 - Change maxrps preset from intrcs + 5 to intrcs + 20
c  mod May 2006 - Change MAXTRS preset to 100
c
C
      PARAMETER (NPARS=7)                                               /* THE NUMBER OF USER PARAMETERS
      DIMENSION SCR(111),LSCR(111),ISCR(111)
      INTEGER*2 ISCR
      INTEGER FIS,SINC,FTR,TRINC,forgat,decimf,order
      COMMON /INPUT/ IPARUN,NLISTS
      CHARACTER*6 NAMES(NPARS)
      CHARACTER*1 TYPES(NPARS)
      DIMENSION LENGTH(NPARS)
      CHARACTER*80 TOKEN
      DIMENSION VALS(NPARS),LVALS(NPARS)
      EQUIVALENCE (VALS(1),LVALS(1))
      COMMON /EDITS/ IERROR,IWARN,IRUN,NOW,ICOMPT
      COMMON /READT/ILUN,NUMHDR,NUMDAT,I1,I2,INTRCS
      COMMON /GEOM/ IGUNIT
      COMMON /TCOL/ LFRP,LRPINC,MWRDS,IOUNIT,NAXRPS,NAXTRS,NINTRS
      INTEGER FRP,RPINC
C
C
      EQUIVALENCE (FRP,LVALS(1)),
     2            (RPINC,LVALS(2)),
     3            (NWRDS,LVALS(3)),
     4            (MAXRPS,LVALS(4)),
     5            (MAXTRS,LVALS(5)),
     6            (MINTRS,LVALS(6)),
     7            (LPRINT,LVALS(7))
      DATA NAMES/'FRP   ','RPINC ','NWRDS ','MAXRPS',
     *           'MAXTRS','MINTRS','LPRINT'/
      DATA LENGTH/3,5,5,4*6/
      DATA TYPES/7*'L'/
C
C
C       SET THE PRESETS
C
      FRP=32767
      RPINC=1
      NWRDS = 0
      MAXRPS=INTRCS+20
      MAXTRS = 100
      MINTRS=1
      NLISTS=0
      NS=0
C****
C****   THE CURRENT COMMAND LINE IN THE SYSTEM BUFFER MAY HAVE THE PARAMETERS.
C****   GET A PARAMETER LIST FROM THE USER.
C****
      NTOKES=1
  100 CONTINUE
      CALL GETOKE(TOKEN,NCHARS)                                          /* GET A TOKEN FROM THE USER PARAMETER LINE
      CALL UPCASE(TOKEN,NCHARS)                                          /* CONVERT THE TOKEN TO UPPERCASE
      IF(NCHARS.GT.0) GO TO 150
      IF(NOW.EQ.1) PRINT 140
  140 FORMAT(' <  ENTER PARAMETERS  >')
      CALL RDLINE                                                        /* GET ANOTHER USER PARAMETER LINE
      NTOKES=0
      GO TO 100
  150 CONTINUE
      NTOKES=NTOKES+1
      DO 190 I=1,NPARS                                                  /* SEE IF IT IS A PARAMETER NAME
      LEN=LENGTH(I)                                                      /* GET THE LEGAL PARAMETER NAME LENGTH
      IPARAM=I                                                          /* SAVE THE INDEX
      IF(TOKEN(1:NCHARS).EQ.NAMES(I)(1:LEN).AND.NCHARS.EQ.LEN) GO TO 200
  190 CONTINUE                                                          /* STILL LOOKING FOR THE NAME
      IF(TOKEN(1:NCHARS).EQ.'END'.AND.NCHARS.EQ.3) GO TO 1000            /* END OF PARAM LIST?
      IF(NS.NE.0) GO TO 230
      PRINT 191, TOKEN(1:NCHARS)
  191 FORMAT(' ***  ERROR  *** GATHER DOES NOT HAVE A PARAMETER ',
     *  'NAMED ',A10)
      IERROR=IERROR+1
      GO TO 100
C****
C****    FOUND THE PARAMETER NAME, NOW FIND THE VALUE
C****
  200 CONTINUE
      NPARAM=IPARAM
  210 CONTINUE                                                           /*  NOW FIND THE VALUE
      CALL GETOKE(TOKEN,NCHARS)
      CALL UPCASE(TOKEN,NCHARS)
      NTOKES=NTOKES+1
      IF(NCHARS.GT.0) GO TO 230                                         /* END OF LINE?
      IF(NOW.EQ.1) PRINT 140                                            /* THIS ALLOWS A PARAMETER TO BE ON A DIFFERENT LINE FROM THE NAME
      CALL RDLINE                                                        /* GET ANOTHER LINE
      NTOKES=0
      GO TO 210
  230 CONTINUE
      IF(TYPES(NPARAM).NE.'A') GO TO 240
      IF(NAMES(NPARAM).EQ.'ADDWB'.AND.TOKEN(1:NCHARS).EQ.'YES')
     *    IADDWB=1
      GO TO 100
  240 CONTINUE
      CALL DCODE(TOKEN,NCHARS,AREAL,ISTAT)                              /* TRY AND DECODE IT
      IF(ISTAT.EQ.2) GO TO 420                                          /* =2 MEANS IT IS A NUMERIC
      IERROR=IERROR+1                                                    /* DCODE PRINTED AN ERROR
      GO TO 100
  420 IF(TYPES(NPARAM).EQ.'L') GO TO 500
      GO TO 100
  500 CONTINUE                                                          /*  32 BIT INTEGER VALUES
      LVALS(NPARAM)=AREAL
      GO TO 100
C****
C****   FINISHED A LIST, NOW DO THE ERROR AND VALIDITY CHECKS
C****
 1000 CONTINUE
      LFRP=FRP                                                           /* FORTRAN DOESN'T ALLOW EQUIVALENCING TO COMMON!!
      LRPINC=RPINC
      MWRDS=NWRDS
      NAXRPS=MAXRPS
      NAXTRS=MAXTRS
      NINTRS=MINTRS
C****
      NS=0
 2020 CALL GETOKE(TOKEN,NCHARS)                                         ! GET THE NEXT TOKEN
      CALL UPCASE(TOKEN,NCHARS)
      NTOKES=NTOKES+1
      IF(NCHARS.GT.0) GO TO 2030                                        ! WAS IT THE END OF A LINE?
      IF(NOW.EQ.1) PRINT 140
      CALL RDLINE                                                       ! GET ANOTHER LINE
      NTOKES=0
      GO TO 2020
 2030 IF(TOKEN(1:NCHARS).NE.'END'.OR.NCHARS.NE.3) GO TO 150
      RETURN                                                            !  FINISHED ALL OF THE PARAMETERS!!!
      END

max_stars_repo_name: apthorpe/sioseis

max_issues_repo_name: henkart/sioseis

max_forks_repo_name: henkart/sioseis


text:
! F2PY experiment
! http://stackoverflow.com/a/12200671/597371

subroutine minmaxr(a,n,amin,amax)
    implicit none
    !f2py intent(hidden) :: n
    !f2py intent(out) :: amin,amax
    !f2py intent(in) :: a
    integer n, i
    real a(n),amin,amax,acurr
    real :: x = 0

    if (n > 0) then
        amin = a(1)
        amax = a(1)
        do i=2, n
            acurr = a(i)
            if (acurr > amax) then
                amax = acurr
            elseif (acurr < amin) then
                amin = acurr
            endif
        enddo
    else
        ! set return values to (inf, -inf)
        amin = -log(x)
        amax = -amin
    endif
end subroutine minmaxr

! Calculate Expected Mutual Information given margins of RxC table
!
! For the sake of numeric precision, the resulting value is *not* normalized
! by N.
!
! License: BSD 3 clause
!
! .. codeauthor:: Robert Layton <robertlayton@gmail.com>
! .. codeauthor:: Corey Lynch <coreylynch9@gmail.com>
! .. codeauthor:: Eugene Scherba <escherba@gmail.com>
!
! List of changes (Eugene Scherba, 10/2/2015):
!
! 1) Removed/rewritten the lines that were creating RxC intermediate NumPy
! arrays which resulted in the O(n^2) memory requirement. Instead, the
! intermediate values are now calculated inside the loop, which may be
! slightly less efficient for small data sizes, but has huge advantages for
! large or even moderately sized data. This change reduces the memory
! requirements of this code from O(n^2) to O(n).
!
! 2) Removed normalization by N from the calculation. It is actually not
! needed to normalize by N if we also don't normalize the input MI value
! (in the calculation of the adjusted score which is ``MI - E(MI) / MI_max
! - E(MI)``, the N value cancels out). Not normalizing the EMI calculation
! by N avoids having to perform lots of tiny floating point increments to
! the EMI aggregate value and thus improves numeric accuracy, especially
! for small values of EMI.
!
! 3) A Fortran 90 version
!
subroutine emi_from_margins(a,R,b,C,emi)
    implicit none
    !f2py intent(hidden) :: R, C
    !f2py intent(out) :: emi
    !f2py intent(in) :: a, b
    integer(kind=8) R, C, i, j, nij, N, N1, N3, max_ab, ai_1, bj_1, N3_ai_1, N3_ai_bj_1

    real(kind=8) emi, log_ai, log_ab_outer_ij, outer_sum, gln_ai_Nai_Ni
    real(kind=8) log_a(R), log_b(C), gln_ai_Nai_N(R), gln_b_Nb(C)

    integer(kind=8) a(R), b(C), a1(R), b1(C)

    real(kind=8), dimension(:), allocatable :: nijs, gln_nij, log_Nnij
    real(kind=8) :: x = -1.0

    log_a = dlog(dble(a))
    log_b = dlog(dble(b))

    N = sum(a)
    if (N /= sum(b)) then
        ! set return value to nan
        emi = dsqrt(x)
        return
    endif

    ! There are three major terms to the EMI equation, which are multiplied to
    ! and then summed over varying nij values.

    ! term1 is nijs.
    ! While nijs[0] will never be used, having it simplifies the indexing.
    max_ab = max(maxval(a), maxval(b))

    allocate(nijs(max_ab + 1))
    do nij=1, max_ab + 1
        nijs(nij) = nij - 1
    enddo

    nijs(1) = 1.0  ! Stops divide by zero warnings. As its not used, no issue.

    ! term2 is log((N*nij) / (a a b)) == log(N * nij) - log(a * b)
    ! term2 uses log(N * nij)
    allocate(log_Nnij(max_ab + 1))
    log_Nnij = dlog(dble(N)) + dlog(nijs)

    ! term3 is large, and involved many factorials. Calculate these in log
    ! space to stop overflows.
    N1 = N + 1
    N3 = N + 3

    a1 = a + 1
    b1 = b + 1
    gln_ai_Nai_N = dlgama(dble(a1)) + dlgama(dble(N1 - a)) - dlgama(dble(N1))
    gln_b_Nb = dlgama(dble(b1)) + dlgama(dble(N1 - b))

    allocate(gln_nij(max_ab + 1))
    gln_nij = dlgama(nijs + 1.0)

    ! emi itself is a summation over the various values.
    emi = 0.0
    do i=1, R
        ai_1 = a1(i)
        log_ai = log_a(i)
        gln_ai_Nai_Ni = gln_ai_Nai_N(i)
        N3_ai_1 = N3 - ai_1
        do j=1, C
            bj_1 = b1(j)
            log_ab_outer_ij = log_ai + log_b(j)
            outer_sum = gln_ai_Nai_Ni + gln_b_Nb(j)
            N3_ai_bj_1 = N3_ai_1 - bj_1

            do nij=max(1, 1 - N3_ai_bj_1), min(ai_1, bj_1) - 1
                ! Numerators are positive, denominators are negative.
                emi = emi + ( &
                    & dble(nij) * &
                    & (log_Nnij(nij + 1) - log_ab_outer_ij) * &
                    & dexp(outer_sum &
                    &    - gln_nij(nij + 1) &
                    &    - dlgama(dble(ai_1 - nij)) &
                    &    - dlgama(dble(bj_1 - nij)) &
                    &    - dlgama(dble(nij + N3_ai_bj_1)) &
                    & ))
            enddo
        enddo
    enddo

    deallocate(gln_nij)
    deallocate(log_Nnij)
    deallocate(nijs)
end subroutine emi_from_margins

max_stars_repo_name: pombredanne/lsh-hdc

max_issues_repo_name: pombredanne/lsh-hdc

max_forks_repo_name: pombredanne/lsh-hdc


text:
C  --------------------------------------------
C  Do matrix inversion for square matrix A(N,N)
C  Output will be written to A(N,N)
C  --------------------------------------------
      subroutine matrix_inverse(N,A)
      IMPLICIT real*8(a-h,o-z)
      integer*4,intent(in) :: N
      complex*16,intent(in) :: A(N,N)
      complex*16 :: WORK(N*N)
      INTEGER :: LWORK, LDA
      INTEGER :: IPIV(N)

C      INTEGER DeAllocateStatus
C      CHARACTER :: A_filename*20
C      external ZGETRF
C      external ZGETRI

      LDA = N
      LWORK = N*N
C      ALLOCATE (A(LDA,N))
C      ALLOCATE (Ainv(LDA,N))
C      ALLOCATE (WORK(LWORK))
C      ALLOCATE (IPIV(N))

C     DGETRF computes an LU factorization of a general M-by-N matrix A
C     using partial pivoting with row interchanges.

C     Store A in Ainv to prevent it from being overwritten by LAPACK
C     Ainv = A
C      print *,'Before inversion'
C      do i=1,N
C        print *,(A(i,j),j=1,N)
C      enddo

      CALL ZGETRF( N, N, A, LDA, IPIV, INFO )

C      IF(INFO.EQ.0)THEN
C         PRINT '(" LU decomposition successful ")'
C      ENDIF
      IF(INFO.LT.0)THEN
         PRINT '(" LU decomposition:  illegal value ")'
         STOP
      ENDIF
      IF(INFO.GT.0)THEN
         WRITE(*,35)INFO,INFO
35       FORMAT( 'LU decomposition: U(',I4,',',I4,') = 0 ')
      ENDIF

C  DGETRI computes the inverse of a matrix using the LU factorization
C  computed by DGETRF.

      CALL ZGETRI(N, A, N, IPIV, WORK, LWORK, INFO)
C      PRINT '(" ")'
      IF (info.NE.0) THEN
         stop 'Matrix inversion failed!'
      ENDIF

C      DEALLOCATE (A, STAT = DeAllocateStatus)
C      DEALLOCATE (Ainv, STAT = DeAllocateStatus)
C      DEALLOCATE (IPIV, STAT = DeAllocateStatus)
C      DEALLOCATE (WORK, STAT = DeAllocateStatus)

      return
      END


text:
!
! test_ClearSky
!
! Test program for the CRTM Adjoint function for clear sky conditions.
!
!

PROGRAM test_ClearSky

  ! ============================================================================
  ! **** ENVIRONMENT SETUP FOR RTM USAGE ****
  !
  ! Module usage
  USE CRTM_Module
  ! Disable all implicit typing
  IMPLICIT NONE
  ! ============================================================================

  ! ----------
  ! Parameters
  ! ----------
  CHARACTER(*), PARAMETER :: PROGRAM_NAME   = 'test_ClearSky'
  CHARACTER(*), PARAMETER :: COEFFICIENTS_PATH = './testinput/'
  CHARACTER(*), PARAMETER :: RESULTS_PATH = './results/adjoint/'

  ! ============================================================================
  ! 0. **** SOME SET UP PARAMETERS FOR THIS TEST ****
  !
  ! Profile dimensions...
  INTEGER, PARAMETER :: N_PROFILES  = 2
  INTEGER, PARAMETER :: N_LAYERS    = 92
  INTEGER, PARAMETER :: N_ABSORBERS = 2
  INTEGER, PARAMETER :: N_CLOUDS    = 0
  INTEGER, PARAMETER :: N_AEROSOLS  = 0
  ! ...but only ONE Sensor at a time
  INTEGER, PARAMETER :: N_SENSORS = 1

  ! Test GeometryInfo angles. The test scan angle is based
  ! on the default Re (earth radius) and h (satellite height)
  REAL(fp), PARAMETER :: ZENITH_ANGLE = 30.0_fp
  REAL(fp), PARAMETER :: SCAN_ANGLE   = 26.37293341421_fp
  REAL(fp), PARAMETER :: SOURCE_ZENITH_ANGLE = 0.0_fp
  REAL(fp), PARAMETER :: AZIMUTH_ANGLE = 135.0_fp
  ! ============================================================================

  ! ---------
  ! Variables
  ! ---------
  CHARACTER(256) :: Message
  CHARACTER(256) :: Version
  CHARACTER(256) :: Sensor_Id
  INTEGER :: Error_Status
  INTEGER :: Allocate_Status
  INTEGER :: n_Channels
  INTEGER :: l, m
  ! Declarations for Adjoint comparisons
  INTEGER :: n_la, n_ma
  INTEGER :: n_ls, n_ms
  CHARACTER(256) :: atmad_file, sfcad_file
  TYPE(CRTM_Atmosphere_type), ALLOCATABLE :: atm_AD(:)
  TYPE(CRTM_Surface_type)   , ALLOCATABLE :: sfc_AD(:)

  ! ============================================================================
  ! 1. **** DEFINE THE CRTM INTERFACE STRUCTURES ****
  !
  TYPE(CRTM_ChannelInfo_type)             :: ChannelInfo(N_SENSORS)
  TYPE(CRTM_Geometry_type)                :: Geometry(N_PROFILES)
  TYPE(CRTM_Options_type)                 :: Options(N_PROFILES)

  ! Define the FORWARD variables
  TYPE(CRTM_Atmosphere_type)              :: Atm(N_PROFILES)
  TYPE(CRTM_Surface_type)                 :: Sfc(N_PROFILES)
  TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution(:,:)

  ! Define the ADJOINT variables
  TYPE(CRTM_Atmosphere_type)              :: Atmosphere_AD(N_PROFILES)
  TYPE(CRTM_Surface_type)                 :: Surface_AD(N_PROFILES)
  TYPE(CRTM_RTSolution_type), ALLOCATABLE :: RTSolution_AD(:,:)
  ! ============================================================================

  !First, make sure the right number of inputs have been provided
  IF(COMMAND_ARGUMENT_COUNT().NE.1)THEN
     WRITE(*,*) TRIM(PROGRAM_NAME)//': ERROR, ONLY one command-line argument required, returning'
     STOP 1
  ENDIF
  CALL GET_COMMAND_ARGUMENT(1,Sensor_Id)   !read in the value

  ! Program header
  ! --------------
  CALL CRTM_Version( Version )
  CALL Program_Message( PROGRAM_NAME, &
    'Test program for the CRTM Adjoint function for clear sky conditions.', &
    'CRTM Version: '//TRIM(Version) )

  ! Get sensor id from user
  ! -----------------------
  Sensor_Id = ADJUSTL(Sensor_Id)
  WRITE( *,'(//5x,"Running CRTM for ",a," sensor...")' ) TRIM(Sensor_Id)

  ! ============================================================================
  ! 2. **** INITIALIZE THE CRTM ****
  !
  ! 2a. Initialise for the requested sensor
  ! ---------------------------------------
  WRITE( *,'(/5x,"Initializing the CRTM...")' )
  Error_Status = CRTM_Init( (/Sensor_Id/), &
                            ChannelInfo  , &
                            File_Path=COEFFICIENTS_PATH)
  IF ( Error_Status /= SUCCESS ) THEN
    Message = 'Error initializing CRTM'
    CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    STOP 1
  END IF

  ! 2b. Determine the total number of channels
  !     for which the CRTM was initialized
  ! ------------------------------------------
  n_Channels = SUM(CRTM_ChannelInfo_n_Channels(ChannelInfo))
  ! ============================================================================

  ! ============================================================================
  ! 3. **** ALLOCATE STRUCTURE ARRAYS ****
  !
  ! 3a. Allocate the ARRAYS
  ! -----------------------
  ALLOCATE( RTSolution( n_Channels, N_PROFILES ), &
            RTSolution_AD( n_Channels, N_PROFILES ), &
            STAT = Allocate_Status )
  IF ( Allocate_Status /= 0 ) THEN
    Message = 'Error allocating structure arrays'
    CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    STOP 1
  END IF

  ! 3b. Allocate the STRUCTURES
  ! ---------------------------
  ! The input FORWARD structure
  CALL CRTM_Atmosphere_Create( Atm, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS )
  IF ( ANY(.NOT. CRTM_Atmosphere_Associated(Atm)) ) THEN
    Message = 'Error allocating CRTM Atmosphere structure'
    CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    STOP 1
  END IF

  ! The output ADJOINT structure
  CALL CRTM_Atmosphere_Create( Atmosphere_AD, N_LAYERS, N_ABSORBERS, N_CLOUDS, N_AEROSOLS )
  IF ( ANY(.NOT. CRTM_Atmosphere_Associated(Atmosphere_AD)) ) THEN
    Message = 'Error allocating CRTM Atmosphere_AD structure'
    CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    STOP 1
  END IF
  ! ============================================================================

  ! ============================================================================
  ! 4. **** ASSIGN INPUT DATA ****
  !
  ! 4a. Atmosphere and Surface input
  ! --------------------------------
  CALL Load_Atm_Data()
  CALL Load_Sfc_Data()

  ! 4b. Geometry input
  ! ------------------
  ! All profiles are given the same value
  !  The Sensor_Scan_Angle is optional.
  CALL CRTM_Geometry_SetValue( Geometry, &
                               Sensor_Zenith_Angle  = ZENITH_ANGLE, &
                               Sensor_Scan_Angle    = SCAN_ANGLE, &
                               Sensor_Azimuth_Angle = AZIMUTH_ANGLE, &
                               Source_Zenith_Angle  = SOURCE_ZENITH_ANGLE )
  ! ============================================================================

  ! ============================================================================
  ! 5. **** INITIALIZE THE ADJOINT ARGUMENTS ****
  !
  ! 5a. Zero the Adjoint OUTPUT structures
  ! ---------------------------------------
  CALL CRTM_Atmosphere_Zero( Atmosphere_AD )
  CALL CRTM_Surface_Zero( Surface_AD )

  ! 5b. Inintialize the Adjoint INPUT so
  !     that the IR/MW results are dTb/dx
  !     and the visible results are dR/dx
  ! -------------------------------------
  DO l = 1, n_Channels
    IF ( ChannelInfo(1)%Sensor_Type == INFRARED_SENSOR .OR. &
         ChannelInfo(1)%Sensor_Type == MICROWAVE_SENSOR ) THEN
      RTSolution_AD(l,:)%Radiance               = ZERO
      RTSolution_AD(l,:)%Brightness_Temperature = ONE
    ELSE
      RTSolution_AD(l,:)%Radiance               = ONE
      RTSolution_AD(l,:)%Brightness_Temperature = ZERO
    END IF
  END DO
  ! ============================================================================

  ! ============================================================================
  ! 6. **** CALL THE CRTM FORWARD MODEL ****
  !
  Error_Status = CRTM_Adjoint( Atm          , &
                               Sfc          , &
                               RTSolution_AD, &
                               Geometry     , &
                               ChannelInfo  , &
                               Atmosphere_AD, &
                               Surface_AD   , &
                               RTSolution   , &
                               Options = Options )
  IF ( Error_Status /= SUCCESS ) THEN
    Message = 'Error in CRTM Adjoint Model'
    CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    STOP 1
  END IF
  ! ============================================================================

  ! ============================================================================
  ! 7. **** OUTPUT THE RESULTS TO SCREEN ****
  !
  DO m = 1, N_PROFILES
    WRITE( *,'(//7x,"Profile ",i0," output for ",a )') m, TRIM(Sensor_Id)
    DO l = 1, n_Channels
      WRITE( *, '(/5x,"Channel ",i0," results")') RTSolution(l,m)%Sensor_Channel
      ! FWD output
      WRITE( *, '(/3x,"FORWARD OUTPUT")')
      CALL CRTM_RTSolution_Inspect(RTSolution(l,m))
    END DO
    ! ADJOINT output
    WRITE( *, '(/3x,"ADJOINT OUTPUT")')
    CALL CRTM_Surface_Inspect(Surface_AD(m))
    CALL CRTM_Atmosphere_Inspect(Atmosphere_AD(m))
  END DO
  ! ============================================================================

  ! ============================================================================
  ! 9. **** COMPARE Atmosphere_AD and Surface_AD RESULTS TO SAVED VALUES ****
  !
  WRITE( *, '( /5x, "Comparing calculated results with saved ones..." )' )

  ! 9a. Create the output files if they do not exist
  ! ------------------------------------------------
  ! 9a.1 Atmosphere file
  ! ...Generate filename
  atmad_file = RESULTS_PATH//TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Atmosphere.bin'
  ! ...Check if the file exists
  IF ( .NOT. File_Exists(atmad_file) ) THEN
    Message = 'Atmosphere_AD save file does not exist. Creating...'
    CALL Display_Message( PROGRAM_NAME, Message, INFORMATION )
    ! ...File not found, so write Atmosphere_AD structure to file
    Error_Status = CRTM_Atmosphere_WriteFile( atmad_file, Atmosphere_AD, Quiet=.TRUE. )
    IF ( Error_Status /= SUCCESS ) THEN
      Message = 'Error creating Atmosphere_AD save file'
      CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
      STOP 1
    END IF
  END IF
  ! 9a.2 Surface file
  ! ...Generate filename
  sfcad_file = RESULTS_PATH//TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Surface.bin'
  ! ...Check if the file exists
  IF ( .NOT. File_Exists(sfcad_file) ) THEN
    Message = 'Surface_AD save file does not exist. Creating...'
    CALL Display_Message( PROGRAM_NAME, Message, INFORMATION )
    ! ...File not found, so write Surface_AD structure to file
    Error_Status = CRTM_Surface_WriteFile( sfcad_file, Surface_AD, Quiet=.TRUE. )
    IF ( Error_Status /= SUCCESS ) THEN
      Message = 'Error creating Surface_AD save file'
      CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
      STOP 1
    END IF
  END IF

  ! 9b. Inquire the saved files
  ! ---------------------------
  ! 9b.1 Atmosphere file
  Error_Status = CRTM_Atmosphere_InquireFile( atmad_file, &
                                              n_Channels = n_la, &
                                              n_Profiles = n_ma )
  IF ( Error_Status /= SUCCESS ) THEN
    Message = 'Error inquiring Atmosphere_AD save file'
    CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    STOP 1
  END IF
  ! 9b.2 Surface file
  Error_Status = CRTM_Surface_InquireFile( sfcad_file, &
                                           n_Channels = n_ls, &
                                           n_Profiles = n_ms )
  IF ( Error_Status /= SUCCESS ) THEN
    Message = 'Error inquiring Surface_AD save file'
    CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    STOP 1
  END IF

  ! 9c. Compare the dimensions
  ! --------------------------
  IF ( n_la /= 0 .OR. n_ma /= N_PROFILES .OR. &
       n_ls /= 0 .OR. n_ms /= N_PROFILES      ) THEN
    Message = 'Dimensions of saved data different from that calculated!'
    CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    STOP 1
  END IF

  ! 9d. Read the saved data
  ! -----------------------
  ! 9d.1 Atmosphere file
  Error_Status = CRTM_Atmosphere_ReadFile( atmad_file, atm_AD, Quiet=.TRUE. )
  IF ( Error_Status /= SUCCESS ) THEN
    Message = 'Error reading Atmosphere_AD save file'
    CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    STOP 1
  END IF
  ! 9d.2 Surface file
  Error_Status = CRTM_Surface_ReadFile( sfcad_file, sfc_AD, Quiet=.TRUE. )
  IF ( Error_Status /= SUCCESS ) THEN
    Message = 'Error reading Surface_AD save file'
    CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    STOP 1
  END IF

  ! 9e. Compare some Adjoints
  ! --------------------------
  ! 9e.1 Atmosphere
  IF ( ALL(CRTM_Atmosphere_Compare(Atmosphere_AD, atm_AD, n_SigFig=3)) ) THEN
    Message = 'Atmosphere_AD Adjoints are the same!'
    CALL Display_Message( PROGRAM_NAME, Message, INFORMATION )
  ELSE
    Message = 'Atmosphere_AD Adjoints are different!'
    CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    ! Write the current Atmosphere_AD results to file
    atmad_file = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Atmosphere.bin'
    Error_Status = CRTM_Atmosphere_WriteFile( atmad_file, Atmosphere_AD, Quiet=.TRUE. )
    IF ( Error_Status /= SUCCESS ) THEN
      Message = 'Error creating temporary Atmosphere_AD save file for failed comparison'
      CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    END IF
  END IF
  ! 9e.2 Surface
  IF ( ALL(CRTM_Surface_Compare(Surface_AD, sfc_AD, n_SigFig=5)) ) THEN
    Message = 'Surface_AD Adjoints are the same!'
    CALL Display_Message( PROGRAM_NAME, Message, INFORMATION )
  ELSE
    Message = 'Surface_AD Adjoints are different!'
    CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    ! Write the current Surface_AD results to file
    sfcad_file = TRIM(PROGRAM_NAME)//'_'//TRIM(Sensor_Id)//'.Surface.bin'
    Error_Status = CRTM_Surface_WriteFile( sfcad_file, Surface_AD, Quiet=.TRUE. )
    IF ( Error_Status /= SUCCESS ) THEN
      Message = 'Error creating temporary Surface_AD save file for failed comparison'
      CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    END IF
  END IF
  ! ============================================================================

  ! ============================================================================
  ! 8. **** DESTROY THE CRTM ****
  !
  WRITE( *, '( /5x, "Destroying the CRTM..." )' )
  Error_Status = CRTM_Destroy( ChannelInfo )
  IF ( Error_Status /= SUCCESS ) THEN
    Message = 'Error destroying CRTM'
    CALL Display_Message( PROGRAM_NAME, Message, FAILURE )
    STOP 1
  END IF
  ! ============================================================================

  ! ============================================================================
  ! 10. **** CLEAN UP ****
  !
  ! 10a. Deallocate the structures
  ! ------------------------------
  CALL CRTM_Atmosphere_Destroy(atm_AD)
  CALL CRTM_Atmosphere_Destroy(Atmosphere_AD)
  CALL CRTM_Atmosphere_Destroy(Atm)

  ! 10b. Deallocate the arrays
  ! --------------------------
  DEALLOCATE(RTSolution, RTSolution_AD, &
             STAT = Allocate_Status)
  ! ============================================================================

  ! Signal the completion of the program. It is not a necessary step for running CRTM.

CONTAINS

  INCLUDE 'Load_Atm_Data.inc'
  INCLUDE 'Load_Sfc_Data.inc'

END PROGRAM test_ClearSky

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 04.02.05.03/04                        *
C  *    TEST TITLE : Validity of predefined edge           *
C  *                 bundles                               *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

C off/on switch for edge flag and error handling mode
      INTEGER    POFF,     PON
      PARAMETER (POFF = 0, PON = 1)
C type of returned value
      INTEGER    PSET,     PREALI
      PARAMETER (PSET = 0, PREALI = 1)

C parameters for <inquire workstation connection and type>
      INTEGER   SPECWT, SPECON

      INTEGER   IDUM1, IDUM2, IDUM3
      INTEGER   MPCI, NPE, EDI, EDI2, ERBUN
      INTEGER   EDCI, EDFL, EDTY, REDTY, PDEF, PDET, PDCI, NAET
      INTEGER   EDCI2, EDFL2, EDTY2
      INTEGER   EREDCI, ERFLG, EREDTY, SVPDET(5), SVPDCI(5), SVPDEF(5)

      REAL      RDUM1, RDUM2, RDUM3, RDUM4(9)
      REAL      EWSF, EWSF2, GENEW, EREW, NOMEW,MINEW,MAXEW
      REAL      SVPDEW(5), PDEW

      LOGICAL   FLGERR, EDGERR, EWERR, COLERR, RLZERR

      CHARACTER MSG*300

      CALL INITGL ('04.02.05.03/04')

C open PHIGS
      CALL XPOPPH (ERRFIL, MEMUN)
C open workstation
      CALL POPWK  (WKID, CONID, WTYPE)
      CALL PQWKC  (WKID, ERRIND, SPECON, SPECWT)
      CALL CHKINQ ('pqwkc', ERRIND)

C  mpci = maximum predefined color index
      CALL PQCF (SPECWT, ERRIND, IDUM1, IDUM2, MPCI, RDUM4)
      CALL CHKINQ ('pqcf', ERRIND)
C  Since color indices are numbered from 0, the maximum index is
C  one less than the number of predefined indices.
      MPCI = MPCI - 1
C  npe = number of predefined entries
C  naet = number of available edgetypes
C  nomew = nominal edgewidth
C  minew = minimum edgewidth
C  maxew = maximum edgewidth
      CALL PQEDF (SPECWT, 0, ERRIND, NAET, IDUM1, IDUM2,
     1            NOMEW, MINEW, MAXEW, NPE)
      CALL CHKINQ ('pqedf', ERRIND)

C  edge flag error found yet
      FLGERR = .FALSE.
C  edgetype error found yet
      EDGERR = .FALSE.
C  edgewidth error found yet
      EWERR  = .FALSE.
C  color error found yet
      COLERR = .FALSE.
C  rlzerr = realized error found yet
      RLZERR = .FALSE.

      CALL SETMSG ('19 20 23', '<Inquire predefined edge '        //
     1             'representation> should be able to access as ' //
     2             'many entries as reported by <inquire edge '   //
     3             'facilities>.')

      DO 100 EDI = 1, NPE
         CALL PQPEDR (SPECWT, EDI, ERRIND, EDFL, EDTY, EWSF, EDCI)
         IF (ERRIND .NE. 0) THEN
            CALL FAIL
            GOTO 150
         ENDIF

         IF (EDFL .LT. POFF .OR. EDFL .GT. PON) THEN
            FLGERR = .TRUE.
            ERFLG = EDFL
         ENDIF

         DO 50 EDI2 = 1, NAET
            CALL PQEDF (SPECWT, EDI2, ERRIND, IDUM1, REDTY, IDUM2,
     1                  RDUM1, RDUM2, RDUM3, IDUM3)
            CALL CHKINQ ('pqedf', ERRIND)
            IF (EDTY .EQ. REDTY) GOTO 60
50       CONTINUE
         EDGERR = .TRUE.
         EREDTY = EDTY

60       CONTINUE
         GENEW = EWSF * NOMEW
         IF (GENEW .LT. MINEW .OR. GENEW .GT. MAXEW) THEN
            EWERR = .TRUE.
            EREW  = EWSF
         ENDIF

         IF (EDCI .LT. 0 .OR. EDCI .GT. MPCI) THEN
            COLERR = .TRUE.
            EREDCI = EDCI
         ENDIF

         CALL PSEDR (WKID, 3, EDFL, EDTY, EWSF, EDCI)
         CALL PQEDR (WKID, 3, PREALI, ERRIND, EDFL2,EDTY2,EWSF2,EDCI2)
         CALL CHKINQ ('pqedr', ERRIND)
         IF (EDFL.EQ.EDFL2 .AND. EDTY.EQ.EDTY2 .AND.
     1       EWSF.EQ.EWSF2 .AND. EDCI.EQ.EDCI2) THEN
C           OK
         ELSE
            RLZERR = .TRUE.
            ERBUN  = EDI
         ENDIF

100   CONTINUE
      CALL PASS

      CALL SETMSG ('20 29', 'The edge flag in predefined edge ' //
     1             'bundles should be either OFF or ON.')
      IF (FLGERR) THEN
         CALL FAIL
         WRITE (MSG, '(A,I5,A)') 'Undefined edge flag', ERFLG,
     1                           ' found in predefined bundle.'
         CALL INMSG (MSG)
      ELSE
         CALL PASS
      ENDIF

      CALL SETMSG ('20 26', 'All edgetypes in predefined edge '   //
     1             'bundles should themselves be in the list of ' //
     2             'available edgetypes.')
      IF (EDGERR) THEN
         CALL FAIL
         WRITE (MSG, '(A,I5,A)') 'Undefined edgetype #', EREDTY,
     1                           ' found in predefined bundle.'
         CALL INMSG (MSG)
      ELSE
         CALL PASS
      ENDIF

      CALL SETMSG ('20 27', 'All edgewidth scale factors in '   //
     1             'predefined edge bundles should generate '   //
     2             'edgewidths within the minimum and maximum ' //
     3             'edgewidths.')
      IF (EWERR) THEN
         CALL FAIL
         WRITE (MSG,'(A,F11.4,A)')
     1          'Unrealizable edgewidth scale factor, ',
     2           EREW, ', found in predefined bundle.'
         CALL INMSG (MSG)
      ELSE
         CALL PASS
      ENDIF

      CALL SETMSG ('20 28', 'All color indices in predefined edge '  //
     1             'bundles should themselves be in the predefined ' //
     2             'color table.')
      IF (COLERR) THEN
         CALL FAIL
         WRITE (MSG, '(A,I5,A)') 'Undefined color #', EREDCI,
     1                           ' found in predefined bundle.'
         CALL INMSG (MSG)
      ELSE
         CALL PASS
      ENDIF

      CALL SETMSG ('15 26 27 28 29', 'All aspects of all predefined ' //
     1             'edge bundles should be reported as realizable '   //
     2             'by <inquire edge representation>.')
      IF (RLZERR) THEN
         CALL FAIL
         WRITE (MSG, '(A,I5,A)') 'Predefined bundle #', ERBUN,
     1                           ' not realized as set.'
         CALL INMSG (MSG)
      ELSE
         CALL PASS
      ENDIF

150   CONTINUE

      WRITE (MSG, '(A,I5,A)') 'PHIGS reported ', NPE,
     1                        ' predefined edge bundles.'
      CALL INMSG (MSG)
      CALL SETMSG ('19 21', '<Inquire edge facilities> should ' //
     1             'report at least 5 predefined entries.')
      IF (NPE .LT. 5) THEN
         CALL FAIL
         GOTO 666
      ELSE
         CALL PASS
      ENDIF

      CALL SETMSG ('20 22', 'The first five predefined edge ' //
     1             'bundle entries should all be mutually '   //
     2             'distinguishable.')

C  pdef = predefined edge flag
C  pdet = predefined edgetype
C  pdew = predefined edgewidth scale factor
C  pdci = predefined edge color index
      DO 200 EDI = 1, 5
C  get attributes for current bundle
         CALL PQPEDR (SPECWT, EDI, ERRIND, PDEF, PDET, PDEW, PDCI)
         CALL CHKINQ ('pqpedr', ERRIND)
C  compare current to all previous bundles
         DO 250 EDI2 = 1, EDI-1
            IF (SVPDEF (EDI2) .EQ. PDEF  .AND.
     1          SVPDET (EDI2) .EQ. PDET  .AND.
     2          SVPDEW (EDI2) .EQ. PDEW  .AND.
     3          SVPDCI (EDI2) .EQ. PDCI)     THEN
               CALL FAIL
               WRITE (MSG, '(A,I2,A,I2,A)') 'Edge bundle #',
     1                EDI, ' is indistinguishable from #', EDI2, '.'
               CALL INMSG (MSG)
               GOTO 666
            ENDIF
250      CONTINUE
C  save current attributes
         SVPDEF (EDI) = PDEF
         SVPDET (EDI) = PDET
         SVPDEW (EDI) = PDEW
         SVPDCI (EDI) = PDCI
200   CONTINUE
      CALL PASS

666   CONTINUE
      CALL ENDIT
      END


text:
great(u1,w1).
great(e1,s1).
great(d1,l1).
great(hh1,p1).
great(d1,ee1).
great(ff1,p1).
great(u1,aa1).
great(z1,s1).
great(y1,cc1).
great(z1,m1).
great(x1,aa1).
great(q1,ee1).
great(ii1,h1).
great(l1,z1).
great(j1,kk1).
great(g1,u1).
great(hh1,s1).
great(a1,r1).
great(ff1,i1).
great(hh1,y1).
great(p1,t1).
great(kk1,s1).
great(aa1,w1).
great(k1,a1).
great(dd1,m1).
great(cc1,o1).
great(ii1,aa1).
great(e1,p1).
great(f1,m1).
great(e1,g1).
great(k1,g1).
great(b1,w1).
great(f1,cc1).
great(ff1,m1).
great(d1,z1).
great(cc1,p1).
great(c1,t1).
great(g1,t1).
great(ff1,cc1).
great(k1,t1).
great(a1,dd1).
great(dd1,s1).
great(a1,g1).
great(a1,aa1).
great(b1,o1).
great(d1,dd1).
great(ii1,p1).
great(kk1,bb1).
great(aa1,o1).
great(a1,z1).
great(q1,ff1).
great(h1,ll1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
      PROGRAM vec2los
C----
C Take three-dimensional N-E-Z displacements and compute line-of-sight
C displacements given the azimuth (from north) and inclination (from
C horizontal) of the view
C----
      IMPLICIT NONE
      REAL*8 pi,d2r
      PARAMETER (pi=4.0d0*atan(1.0d0),d2r=pi/180.0d0)
      REAL*8 uLOS,stlo,stla,stdp,uN,uE,uZ
      REAL*8 az,incl,cosaz,cosinc,sinaz,sininc
      CHARACTER*40 ifile,ofile
      LOGICAL ex

C----
C Parse command line
C----
      call gcmdln(az,incl,ifile,ofile)
      if (az.lt.-1000.0d0.or.incl.lt.-1000.0d0) then
          write(*,*) '!! Error: azimuth or inclination are undefined'
          call usage('!! Use -a AZ or -i INCL')
      endif
      if (ifile.eq.'none') then
          write(*,*) '!! Error: Input displacement vector file ',
     1               'unspecified'
          call usage('!! Use -f IFILE to specify input file')
      else
          inquire(file=ifile,EXIST=ex)
          if (.not.ex) call usage('!! Error: no input file: '//ifile)
      endif
      if (ofile.eq.'none') then
          write(*,*) '!! Error: No output file specified'
          call usage('!! Use -o OFILE to specify LOS output file')
      endif
      cosaz = dcos(az*d2r)
      sinaz = dsin(az*d2r)
      cosinc = dcos(incl*d2r)
      sininc = dsin(incl*d2r)
      !print *,'AZ',az
      !print *,'INCL',incl
      !print *,'IFILE ',ifile
      !print *,'OFILE ',ofile

C----
C Control files
C----
      open (unit=11,file=ifile,status='old')
      open (unit=12,file=ofile,status='unknown')

C----
C Convert NEZ to LOS
C----
   11 read (11,*,end=12) stlo,stla,stdp,uE,uN,uZ
          uLOS = uN*cosaz*cosinc
     1             + uE*sinaz*cosinc
     2                 - uZ*sininc
          write (12,9991) stlo,stla,stdp,uLOS
          goto 11
   12 continue

 9991 format(3F12.4,X,3E12.4)

      END

C======================================================================C

      SUBROUTINE gcmdln(az,incl,ifile,ofile)
      IMPLICIT none
      REAL*8 az,incl
      CHARACTER*40 tag,ifile,ofile
      INTEGER i,narg
      az = -99999.0d0
      incl = -99999.0d0
      ifile = 'none'
      ofile = 'none'
      narg = iargc()
      if (narg.eq.0) then
          call usage('!! Error: no command line arguments specified')
      endif
      i = 0
   11 i = i + 1
      if (i.gt.narg) goto 12
          call getarg(i,tag)
          if (tag(1:2).eq.'-a') then
              i = i + 1
              call getarg(i,tag)
              read (tag,'(BN,F12.0)') az
          elseif (tag(1:2).eq.'-i') then
              i = i + 1
              call getarg(i,tag)
              read (tag,'(BN,F12.0)') incl
          elseif (tag(1:2).eq.'-f') then
              i = i + 1
              call getarg(i,ifile)
          elseif (tag(1:2).eq.'-o') then
              i = i + 1
              call getarg(i,ofile)
          elseif (tag(1:2).eq.'-h') then
              call usage(' ')
          endif
          goto 11
   12 continue
      RETURN
      END

C----------------------------------------------------------------------C

      SUBROUTINE usage(str)
      IMPLICIT none
      INTEGER lstr
      CHARACTER str*(*)
      if (str.ne.' ') then
          lstr = len(str)
          write(*,*) str(1:lstr)
          write(*,*)
      endif
      write(*,*)
     1 'Usage: vec2los -a AZ -i INC -f IFILE -o OFILE [-h]'
      write(*,*)
      write(*,*)
     1 '-a AZ     Look azimuth (CW from N)'
      write(*,*)
     1 '-i INC    Look inclination (from horizontal)'
      write(*,*)
     1 '-f IFILE  Input vector displacement file (stlo stla stdp uE ',
     2                                                'uN uZ)'
      write(*,*)
     1 '-o OFILE  Output LOS displacement file (stlo stla stdp uLOS)'
      write (*,*)
     1 '-h        Online help (this message)'
      write (*,*)
      STOP
      END


text:
      SUBROUTINE ZROT( N, CX, INCX, CY, INCY, C, S )
*
*  -- LAPACK auxiliary routine (version 3.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            INCX, INCY, N
      DOUBLE PRECISION   C
      COMPLEX*16         S
*     ..
*     .. Array Arguments ..
      COMPLEX*16         CX( * ), CY( * )
*     ..
*
*  Purpose
*  =======
*
*  ZROT   applies a plane rotation, where the cos (C) is real and the
*  sin (S) is complex, and the vectors CX and CY are complex.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The number of elements in the vectors CX and CY.
*
*  CX      (input/output) COMPLEX*16 array, dimension (N)
*          On input, the vector X.
*          On output, CX is overwritten with C*X + S*Y.
*
*  INCX    (input) INTEGER
*          The increment between successive values of CY.  INCX <> 0.
*
*  CY      (input/output) COMPLEX*16 array, dimension (N)
*          On input, the vector Y.
*          On output, CY is overwritten with -CONJG(S)*X + C*Y.
*
*  INCY    (input) INTEGER
*          The increment between successive values of CY.  INCX <> 0.
*
*  C       (input) DOUBLE PRECISION
*  S       (input) COMPLEX*16
*          C and S define a rotation
*             [  C          S  ]
*             [ -conjg(S)   C  ]
*          where C*C + S*CONJG(S) = 1.0.
*
* =====================================================================
*
*     .. Local Scalars ..
      INTEGER            I, IX, IY
      COMPLEX*16         STEMP
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          DCONJG
*     ..
*     .. Executable Statements ..
*
      IF( N.LE.0 )
     $   RETURN
      IF( INCX.EQ.1 .AND. INCY.EQ.1 )
     $   GO TO 20
*
*     Code for unequal increments or equal increments not equal to 1
*
      IX = 1
      IY = 1
      IF( INCX.LT.0 )
     $   IX = ( -N+1 )*INCX + 1
      IF( INCY.LT.0 )
     $   IY = ( -N+1 )*INCY + 1
      DO 10 I = 1, N
         STEMP = C*CX( IX ) + S*CY( IY )
         CY( IY ) = C*CY( IY ) - DCONJG( S )*CX( IX )
         CX( IX ) = STEMP
         IX = IX + INCX
         IY = IY + INCY
   10 CONTINUE
      RETURN
*
*     Code for both increments equal to 1
*
   20 CONTINUE
      DO 30 I = 1, N
         STEMP = C*CX( I ) + S*CY( I )
         CY( I ) = C*CY( I ) - DCONJG( S )*CX( I )
         CX( I ) = STEMP
   30 CONTINUE
      RETURN
      END

max_issues_repo_name: GEHC-Surgery/ITK

max_forks_repo_name: GEHC-Surgery/ITK


text:
Classical and Contemporary Dance Classes
with Timothy Nutter

Timothys teaching style emphasizes aspects of dance relating to theater and music: emotion, gesture, tempo and musicality. He also believes in encouraging students to think and evaluate for themselves, to draw connections between dance and other areas of experience, and always to continue making adjustments.

Contact Timothy by email at: timothynutterdance@gmail.com

http://www.youtube.com/watch?vDIwrGaOyRMI Click here for a video of a recent dance and music collaboration with guitarist Matthew Grasso.

Very Beginning Adult Ballet
My apologies, but this class is now cancelled! My best wishes to all the adult dancers in Davis, and good luck with your studies.

Here is a beginning ballet class just for adults. If you have never tried ballet before, this is the class to get you started. We will cover all the fundamental skills, exercises and vocabulary of basic ballet technique, but at a slow pace, ideal for beginning adults. It’s a great way to exercise the whole body, while you enjoy learning about a creative and beautiful dance tradition.

The Time: 	Fridays, 11:00am–12:30pm (Please note the new later time slot will begin on November 8th.)

The Place:	http://www.davisholistic.50webs.com Davis Holistic Health Center, 1403 5th Street, Suite B

First class will be October 18th. For now the fee will be $15 per 90 minute class, paid monthly, cash or personal check. (October fee will be $30, November $75, etc.)

Ballet Partnering

The Place: http://www.davisholistic.50webs.com Davis Holistic Health Center, 1403 5th Street, Suite B

This class is suspended indefinitely! My apologies for the inconvenience.

The Time: Sundays 11am–12:30pm

For now the fee will remain what it was at FDF, $15 per 90 minute class. Please bring cash or check, as I will not be able to take debit or credit cards.

The Davis Holistic Health Center is on the northeastern corner of the intersection of L Street and 5th Street. The building is east of the corner gas station, along 5th Street. To get there by car, turn into the gas station, then drive around the buildings behind the station to get to the parking lot for the Center and the Fitness Garage.

Biography
TimothyNutter Timothy Nutter attended college at SUNY Purchase and Antioch College. At Purchase he studied in the Conservatory of Dance, while producing his own dance shows and running an independent improv dance troupe. During this time he also studied in New York City, taking class with contemporary artists such as Wil Swanson, Keith Thompson and Kathleen Fisher (all formerly of Trisha Brown Dance Company), and learning contact improvisation in the basement of Judson Church.

At Antioch Timothy studied philosophy, linguistics and drama. He graduated in 2004 with a B.A. in philosophy. Timothy attended graduate school in 2005 at Villanova University in Philadelphia, where he studied history of philosophy. In 2012, he again attended graduate school, this time working on the philosophy of biology at Exeter University in England.

Timothy has trained in a number of different dance traditions, but mainly ballet, modern, postmodern, improvisation, and ballroom. Teachers in the Davis and Sacramento area he has worked with include: Hanneke Lohse, Pamela Trokanski, Barbara Nicholas, Gwen Burton, Lisa Applegate, and Cara Rains. Two of Timothys philosophy professors were Scott Warren at Antioch and James Wetzel at Villanova. At Purchase he studied music with the pianist Richard CameronWolfe and dance composition with Neil Greenberg. In New York, he studied acting at the Terry Schreiber Studio and Theatre.

Timothy has also begun a small business as a The Choice of Words writing tutor.

Timothy regularly performs his songs on guitar at the Monticello Monticello Restaurant open mic nights (Wednesdays at 8:30pm).

Follow the truth, wherever it may lead. –Thomas Jefferson