text:
C    %W% %G%
      subroutine wrthis
C
C     This subroutine writes solution data to the history file FOR
C     It is called by SOLN and CNTRL.
C
      include 'tspinc/params.inc'
      include 'tspinc/newton.inc'
      include 'tspinc/buslod.inc'
      include 'tspinc/blkcom1.inc'
      include 'tspinc/cntrl2.inc'
      include 'tspinc/prt.inc'
      include 'tspinc/param.inc'
      include 'tspinc/ectba.inc'
      include 'tspinc/mdctbl.inc'
      include 'tspinc/gamma.inc'
      include 'tspinc/dcmodd.inc'
      include 'tspinc/machd1.inc'
      include 'tspinc/machd2.inc'
      include 'tspinc/outaux.inc'
      include 'tspinc/busvolt.inc'
      include 'tspinc/znox.inc'
      include 'tspinc/znox2.inc'
      include 'tspinc/vy1.inc'
      include 'tspinc/vym.inc'
      include 'tspinc/ecsind.inc'
      include 'tspinc/bname.inc'
      include 'tspinc/buskv.inc'
      include 'tspinc/bcur.inc'
      include 'tspinc/igentn.inc'
      include 'tspinc/spare1.inc'
      include 'tspinc/vymn.inc'

C     Functions

      logical dbghere

C     Local variables

      dimension temps1(390), temps3(900)
      dimension tempn(2, 13)
      character*8 name, name1, name2
      character*1 id
      logical debug, first

      save

      data radeg/57.29578/
      data debug/.true./
      data to_old, e_old, f_old / 3 * 0.0 /

      data first / .true./

C
C     Begin   Begin   Begin   Begin   Begin   Begin   Begin   Begin
C
      if (first) then

        write (98, 94) 'TimeBase=Cycles'
        write (98, 94) '8'
        write (98, 94) 'Time VNET_MAG VNET_ANG VMAC_MAG VMAC_ANG INET_MA
     &G INET_ANG IMAC_MAG'
   94   format (a)

        write (99, 94) 'TimeBase=Cycles'
        write (99, 94) '9'
        write (99, 94) 'Time IMAC_ANG FD FQ IND_FREQ SLIP PELECT PMECH F
     &REQ'
        first = .false.
      endif

      debug = dbghere('WRTHIS  ')
      if (debug) then
        call dbgeko('WRTHIS - at start')
        call dbgwrf('  TO /current time/ = ', to)
      endif
C
C     CALCULATE BUS LOAD BUSP - REAL LOAD IN MW'S
C     BUSQ -  REACTIVE LOAD IN MVAR'S
C
C     *** CHANGED BUS LOAD OPTION TEST FROM ISPSW TO IBLOD

      if (iblod .ne. 0) then

        do itrr = 1, nmx
          name = bname(itrr)
          curr =  - busp(itrr) + eyr(itrr)*yreal(itrr) - eyi(itrr)
     &     *yimag(itrr)
          curi =  - busq(itrr) + eyi(itrr)*yreal(itrr) + eyr(itrr)
     &     *yimag(itrr)
          busp(itrr) = (eyr(itrr)*curr+eyi(itrr)*curi)*bmva
          busq(itrr) = (curr*eyi(itrr)-eyr(itrr)*curi)*bmva
          if (keybrd(0) .ne. 0) then
            write (outbuf, 10000) name, busp(itrr), busq(itrr)
10000       format (5x, a8, 2(2x, f9.2))
            call prtout(1)
          endif
        enddo
      endif
      ifirst = 5
      jfirst = 1 + 6*isg
      if (iblod .ne. 0) jfirst = 1 + 7*isg
      if (ispsw .ne. 0) jfirst = 1 + 8*isg

C     ADDITIONAL CODE FOR BOTH IBLOD,ISPSW .NE. 0

      if (iblod .ne. 0 .and. ispsw .ne. 0) jfirst = 1 + 9*isg
C
C     Write generator data to solution file (using unformatted,
c     direct-access writes).

c     Record No. 130 = time step

      irectp = 130
      idesc = 130
      irecln = 1
      if (debug) then
        call dbgeko2('WRTHIS - saving simulation time to ',
     &   'history file.')
        call dbgwri('  IRECTP /record type/ = ', irectp)
        call dbgwri('  IDESC /rec descrip/  = ', idesc)
        call dbgwri('  IRECLN /rec length/  = ', irecln)
      endif
      call puthisrf(irectp, idesc, irecln, to)

c     Record No. 134 = gen rotor angles

      irectp = 134
      idesc = 134
      irecln = isg
      if (debug) then
        call dbgeko2('WRTHIS - saving gen rotor angles to ',
     &   'history file.')
        call dbgwri('  IRECTP /record type/ = ', irectp)
        call dbgwri('  IDESC /rec descrip/  = ', idesc)
        call dbgwri('  IRECLN /rec length/  = ', irecln)
      endif
      call puthisrf(irectp, idesc, irecln, angl)

c     Record No. 138 = gen frequency devs

      irectp = 138
      idesc = 138
      irecln = isg
      if (debug) then
        call dbgeko2('WRTHIS - saving gen rotor frequency devs to ',
     &   'history file.')
        call dbgwri('  IRECTP /record type/ = ', irectp)
        call dbgwri('  IDESC /rec descrip/  = ', idesc)
        call dbgwri('  IRECLN /rec length/  = ', irecln)
      endif
      call puthisrf(irectp, idesc, irecln, angp2)

c     Record No. 142 = gen field voltages & currents

      irectp = 142
      idesc = 142
      irecln = isg*2
      if (debug) then
        call dbgeko2('WRTHIS - saving field volts & currents to ',
     &   'history file.')
        call dbgwri('  IRECTP /record type/ = ', irectp)
        call dbgwri('  IDESC /rec descrip/  = ', idesc)
        call dbgwri('  IRECLN /rec length/  = ', irecln)
      endif
      call puthisrf(irectp, idesc, irecln, vfldtn)

c     Record No. 146 = gen mech power

      irectp = 146
      idesc = 146
      irecln = isg
      if (debug) then
        call dbgeko2('WRTHIS - saving gen P_mech to ', 'history file.')
        call dbgwri('  IRECTP /record type/ = ', irectp)
        call dbgwri('  IDESC /rec descrip/  = ', idesc)
        call dbgwri('  IRECLN /rec length/  = ', irecln)
      endif
      call puthisrf(irectp, idesc, irecln, govpwr)

c     Record No. 150 = gen elec power

      irectp = 150
      idesc = 150
      irecln = isg
      if (debug) then
        call dbgeko2('WRTHIS - saving gen P_elec to ', 'history file.')
        call dbgwri('  IRECTP /record type/ = ', irectp)
        call dbgwri('  IDESC /rec descrip/  = ', idesc)
        call dbgwri('  IRECLN /rec length/  = ', irecln)
      endif
      call puthisrf(irectp, idesc, irecln, genp)

c     Record No. 10150 = gen flux

      irectp = 10150
      idesc = 10150
      irecln = isg
      if (debug) then
        call dbgeko2('WRTHIS - saving gen flux ', 'history file.')
        call dbgwri('  IRECTP /record type/ = ', irectp)
        call dbgwri('  IDESC /rec descrip/  = ', idesc)
        call dbgwri('  IRECLN /rec length/  = ', irecln)
      endif
      call puthisrf(irectp, idesc, irecln, epq)

c     Record No. 154 = gen terminal current

      irectp = 154
      idesc = 154
      irecln = isg*2
      if (debug) then
        call dbgeko2('WRTHIS - saving gen term current (R & I) to ',
     &   'history file.')
        call dbgwri('  IRECTP /record type/ = ', irectp)
        call dbgwri('  IDESC /rec descrip/  = ', idesc)
        call dbgwri('  IRECLN /rec length/  = ', irecln)
      endif
      call puthisrf(irectp, idesc, irecln, curn)

c     Record No. 158 = gen voltage reg output

      irectp = 158
      idesc = 158
      irecln = isg
      if (debug) then
        call dbgeko2('WRTHIS - saving gen voltage reg output to ',
     &   'history file.')
        call dbgwri('  IRECTP /record type/ = ', irectp)
        call dbgwri('  IDESC /rec descrip/  = ', idesc)
        call dbgwri('  IRECLN /rec length/  = ', irecln)
      endif
      call puthisrf(irectp, idesc, irecln, regout)

c     Record No. 162 = gen PSS output

      irectp = 162
      idesc = 162
      irecln = isg
      if (debug) then
        call dbgeko2('WRTHIS - saving gen PSS output to ',
     &   'history file.')
        call dbgwri('  IRECTP /record type/ = ', irectp)
        call dbgwri('  IDESC /rec descrip/  = ', idesc)
        call dbgwri('  IRECLN /rec length/  = ', irecln)
      endif
      call puthisrf(irectp, idesc, irecln, supout)
      if (iblod .ne. 0) then

c       Record No. 166 = bus loads (real)

        irectp = 166
        idesc = 166
        irecln = nmx
        if (debug) then
          call dbgeko2('WRTHIS - saving bus loads (real) to ',
     &     'history file.')
          call dbgwri('  IRECTP /record type/ = ', irectp)
          call dbgwri('  IDESC /rec descrip/  = ', idesc)
          call dbgwri('  IRECLN /rec length/  = ', irecln)
        endif
        call puthisrf(irectp, idesc, irecln, busp)

c       Record No. 170 = bus loads (reactive)

        irectp = 170
        idesc = 170
        irecln = nmx
        if (debug) then
          call dbgeko2('WRTHIS - saving bus loads (imag) to ',
     &     'history file.')
          call dbgwri('  IRECTP /record type/ = ', irectp)
          call dbgwri('  IDESC /rec descrip/  = ', idesc)
          call dbgwri('  IRECLN /rec length/  = ', irecln)
        endif
        call puthisrf(irectp, idesc, irecln, busq)
      endif
      if (ispsw .ne. 0) then

c       Record No. 258 = spare points

        irectp = 258
        idesc = 258
        irecln = MAXSP
        if (debug) then
          call dbgeko2('WRTHIS - saving spare points to ',
     &     'history file.')
          call dbgwri('  IRECTP /record type/ = ', irectp)
          call dbgwri('  IDESC /rec descrip/  = ', idesc)
          call dbgwri('  IRECLN /rec length/  = ', irecln)
        endif
        call puthisrf(irectp, idesc, irecln, sparpt)
C       WRITE(L8) (SPARPT(I),I=1,MAXSP)
      endif

c     Record No. 174 = bus voltages (real)

      irectp = 174
      idesc = 174
      irecln = nmx
      if (debug) then
        call dbgeko2('WRTHIS - saving bus voltage (real) to ',
     &   'history file.')
        call dbgwri('  IRECTP /record type/ = ', irectp)
        call dbgwri('  IDESC /rec descrip/  = ', idesc)
        call dbgwri('  IRECLN /rec length/  = ', irecln)
      endif
      call puthisrf(irectp, idesc, irecln, eyr)

c     Record No. 178 = bus voltages (imaginary)

      irectp = 178
      idesc = 178
      irecln = nmx
      if (debug) then
        call dbgeko2('WRTHIS - saving bus voltage (imag) to ',
     &   'history file.')
        call dbgwri('  IRECTP /record type/ = ', irectp)
        call dbgwri('  IDESC /rec descrip/  = ', idesc)
        call dbgwri('  IRECLN /rec length/  = ', irecln)
      endif
      call puthisrf(irectp, idesc, irecln, eyi)
C
C     Output variable admittance values (Printer only)
C
      if (iznmax .gt. 0 .and. keybrd(34) .ne. 0) then
        do i = 1, iznmax
          name1 = bname(iznbus(i))
          bkv1 = buskv(iznbus(i))
          name2 = bname(jznbus(i))
          bkv2 = buskv(jznbus(i))
          id = iznpar(i)
          alfa = (vyalfi(i)+vyalfa(i))*radeg
          if (iznosw(i) .ne. 1) then
            write (outbuf, 10010) name1, bkv1, name2, bkv2, id, alfa,
     &       zngij(i), znbij(i), vymmd(i), to
10010       format ('0', 5x, 2(a8, 1x, f4.0, 1x), a1, ' ALPHA = ',
     &       f6.2, ' GIJ = ', f9.4, ' BIJ = ', f9.4, 2x, ' MOD SIG = ',
     &       f9.4, ' AT ', f7.2, ' CYCLES ')
            call prtout(1)
          endif
          if (iznosw(i) .ne. 3) then
            recp2 = 1./(zngc(i)*zngc(i)+znbc(i)*znbc(i))
            rtot = zngc(i)*recp2
            xtot =  - znbc(i)*recp2
            write (outbuf, 10020) name1, bkv1, name2, bkv2, id, rtot,
     &       xtot, to
10020       format ('0', 5x, 2(a8, 1x, f4.0, 1x), a1, 1x, ' ZNO R = ',
     &       f9.4, ' ZNO X = ', f9.4, 'AT ', f7.2, ' CYCLES ')
            call prtout(1)
          endif
        enddo
      endif
C
C     Output dc quantities (Printer only)
C
      if (.not. (ldc .eq. 0 .and. ldc1 .eq. 0)) then
        i784 = 1
        idc1 = 0
        idc2 = 0
        if (ldc .ne. 0) then
          i1 = 0
          do i = 1, ldc
            i2 = (i-1)*idcl + kdc
            call redecs(ptab(1), i2, 252)
            if (.not. (itab(114) .lt. 1 .or. itab(114) .eq. 5)) then
              imod = itab(118)
              if (keybrd(22) .ne. 0) then
                if (itab(114) .eq. 6) then
                  write (outbuf, 10030) to, x1nmod(imod), y1nmod(imod),
     &             sighim(imod)
                  call prtout(1)
10030             format ('0', 5x, 'TIME = ', f7.2, 5x, 'MOD1 INPUT ',
     &             e13.6, 5x, 'MOD2 INPUT ', e13.6, 5x, 'DUAL MOD OUT '
     &             , e13.6)
                else
                  write (outbuf, 10040) to, xoomod(imod), siglom(imod),
     &             sighim(imod)
                  call prtout(1)
                endif
10040           format ('0', 5x, 'TIME = ', f7.2, 5x, 'MOD INPUT ',
     &           e13.6, 5x, 'LOW LEVEL OUTPUT ', e13.6, 5x,
     &           'HIGH LEVEL OUTPUT ', e13.6)
              endif
            endif
C
C           ITEMPSC(I) = RECTIFIER BUS NAME CHAR*8
C           ITEMPSC(I+1) = INVERTER BUS NAME CHAR*8
C           IPTAB(44) = RECTIFIER BASE KV CODE
C           IPTAB(45) = INVERTER BASE KV CODE
C           TAB(95) = COS ALPHA FOR RECTIFIER
C           TAB(87) = IDC FOR RECTIFIER
C           TAB(110) = VDC FOR RECTIFIER
C           TAB(157) = IDC FOR INVERTER
C           TAB(111) = VDC FOR INVERTER
C           ITAB(114) IS THE MODULATION CODE
C           TAB(121) IS THE EXTINCTION ANGLE FOR THE RECTIFIER
C           TAB(122) IS THE EXTINCTION ANGLE FOR THE INVERTER
C           TAB(115) IS V ALPHA FOR THE RECTIFIER
C           TAB(116) IS V ALPHA FOR THE INVERTER
C           TAB(78) IS V ALPHA PRIME FOR THE RECTIFIER
C           TAB(84) IS V ALPHA PRIME FOR THE INVERTER
C
            temps1(18*i-17) = tab(95)
            temps1(18*i-16) = tab(87)
            temps1(18*i-15) = tab(110)
            temps1(18*i-14) = ptab(25)
            temps1(18*i-13) = tab(121)
            temps1(18*i-12) = tab(115)
            temps1(18*i-11) = tab(78)
            temps1(18*i-10) = 0.0
            temps1(18*i-9) = 0.0
            temps1(18*i-8) = tab(96)
            temps1(18*i-7) = tab(157)
            temps1(18*i-6) = tab(111)
            temps1(18*i-5) = ptab(26)
            temps1(18*i-4) = tab(122)
            temps1(18*i-3) = tab(116)
            temps1(18*i-2) = tab(84)
            temps1(18*i-1) = 0.0
            temps1(18*i) = 0.0
C
C           PRINT WARNING IF ANGLE OF OVERLAP EXCEEDS 60 DEGREES
C
            aover = (temps1(18*i-13)-acos(tab(95)))*57.29578
            if (aover .gt. 60.) then
              ibno = itab(35)
              name = bname(ibno)
              bkv = buskv(ibno)
              write (errbuf(1), 10050) name, bkv, aover
10050         format ('0', 2x, 'ANGLE OF OVERLAP AT ', a8, 1x, f6.2,
     &         ' EXCEEDS 60', ' DEGREES.  ANGLE =  ', f6.2)
              call prterr('W', 1)
            endif
            aover = (temps1(18*i-4)-acos(tab(96)))*57.29578
            if (aover .gt. 60.) then
              ibno = itab(37)
              name = bname(ibno)
              bkv = buskv(ibno)
              write (errbuf(1), 10050) name, bkv, aover
              call prterr('W', 1)
            endif
C
C           STORE MODULATION OUTPUT
C           GAMMA MODULATION STORED WITH INVERTER BUS
C
            if (itab(114) .eq. 5) then
              igam = itab(118)
              temps1(18*i-5) = gama(igam)
            endif

C           LOW LEVEL MODULATION

            if (itab(114) .eq. 1 .or. itab(114) .eq. 2)
     &       temps1(18*i-14) = siglo

C           HIGH LEVEL MODULATION

            if (itab(114) .eq. 3 .or. itab(114) .eq. 4)
     &       temps1(18*i-14) = sighi

C           DUAL FREQUENCY MODULATION

            if (itab(114) .eq. 6) temps1(18*i-14) = tab(177)
          enddo
          idc1 = ldc*18
        endif
        if (ldc1 .ne. 0) then
          ldcc = 1
          n999 =  - 999
          jecs1 = kdc1
          do while (.true.)
            call redecs(dcb, jecs1, lbt)

C           INITIALIZE NFLT FOR FLT  DC BKR TRIP CODE IF REQ

            if (nflt .eq. 999) then
              kflt = 0
              call ritecs(kflt, jecs1+7, 1)
              call ritecs(n999, jecs1+13, 1)
            endif
            ndimj = 3*(ndimh-nterm)/2 + 9*nterm
            call redecs(dcj, necsj, ndimj)

C           PROCESS TERMINAL DATA FIRST

            nterm3 = 9*nterm
            do i = 1, nterm3
              temps1(idc1+i) = dcj(i)
            enddo
            idc1 = idc1 + nterm3

C           PROCESS DC BRANCH DATA

            i1 = 0
            nterm3 = nterm3 + 1
            do i = nterm3, ndimj
              i1 = i1 + 1
              temps3(idc2+i1) = dcj(i)
            enddo
            idc2 = idc2 + ndimj - nterm3 + 1
C
C           LOOP THROUGH TERMINALS AND OUTPUT ANY MODULATION SIGNA
C
            if (keybrd(22) .ne. 0) then
              ind = 0
              do while (.true.)
                ind = ind + 1
                if (ind .gt. nterm) goto 100
                ind1 = lbt + ind
                call redecs(dca, idcb(ind1), msizea)
C
C               PRINT WARNING IF ANGLE OF OVERLAP EXCEEDS 60 DEGRE
C
                itrr = 9*(ind-1)
                aover = (dcj(itrr+5)-acos(dcj(itrr+1)))*57.29578
                if (aover .gt. 60.) then
                  ibno = idca(46)
                  name = bname(ibno)
                  bkv = buskv(ibno)
                  write (errbuf(1), 10050) name, bkv, aover
                  call prterr('W', 1)
                endif
                modcod = idca(78)
                if (.not. (modcod .lt. 1 .or. modcod .eq. 5)) then
                  imod = idca(79)
                  if (modcod .eq. 6) then
                    write (outbuf, 10030) to, x1nmod(imod), y1nmod
     &               (imod), sighim(imod)
                    call prtout(1)
                  else
                    write (outbuf, 10040) to, xoomod(imod), siglom
     &               (imod), sighim(imod)
                    call prtout(1)
                  endif
                endif
              enddo
            endif
C
C           TALLY AND TEST END OF MULTDC CKT
C
  100       ldcc = ldcc + 1
            if (ldcc .gt. ldc1) goto 110
            jecs1 = kdcm(ldcc)
          enddo
        endif
C
C       WRITE VARIABLE ADMITTANCE FIRING ANGLE OVER DC DATA
C
  110   if (iznmax .gt. 0) then
          temps1(1) = cos(vyalfa(1)+vyalfi(1))
          temps1(2) = vymmd(1)
        endif

C       Record No. 182 = D-C Bus data

        irectp = 182
        idesc = 182
        irecln = idc1
        if (debug) then
          call dbgeko2('WRTHIS - saving DC bus data to ',
     &     'history file.')
          call dbgwri('  IRECTP /record type/ = ', irectp)
          call dbgwri('  IDESC /rec descrip/  = ', idesc)
          call dbgwri('  IRECLN /rec length/  = ', irecln)
        endif
        call puthisrf(irectp, idesc, irecln, temps1)

C       Record No. 186 = D-C Bus data

        if (idc2 .ne. 0) then
          irectp = 186
          idesc = 186
          irecln = idc2
          if (debug) then
            call dbgeko2('WRTHIS - saving DC branch data to ',
     &       'history file.')
            call dbgwri('  IRECTP /record type/ = ', irectp)
            call dbgwri('  IDESC /rec descrip/  = ', idesc)
            call dbgwri('  IRECLN /rec length/  = ', irecln)
          endif
          call puthisrf(irectp, idesc, irecln, temps3)
C         WRITE(L8)(TEMPS3(I),I=1,IDC2)
        endif
      endif

C     Record No. 190 = TCSC data points

      if (iznmax .ne. 0) then
        irectp = 190
        idesc = 190
        irecln = 2*iznmax
        if (debug) then
          call dbgeko2('WRTHIS - saving TCSC data to ',
     &     'history file.')
          call dbgwri('  IRECTP /record type/ = ', irectp)
          call dbgwri('  IDESC /rec descrip/  = ', idesc)
          call dbgwri('  IRECLN /rec length/  = ', irecln)
        endif
        do i = 1, iznmax
          if (vyalfi(i) .gt. 0.0 .and. iranityp(i) .eq. 3) then
            zbase = buskv(iznbus(i))**2 / bmva
            temps1(2*i-1) = -xtcscge(i) * zbase
            temps1(2*i) = dangl(i) * 57.2957795
          else
            temps1(2*i-1) = 0.0
            temps1(2*i) = 0.0
          endif
        enddo
        call puthisrf(irectp, idesc, irecln, temps1)
      endif

C     Write an end-of-file marker in case this is the last time ste

      call puthisw(-1)

      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
C$Procedure      PXFORM ( Position Transformation Matrix )

      SUBROUTINE PXFORM ( FROM, TO, ET, ROTATE )

C$ Abstract
C
C     Return the matrix that transforms position vectors from one
C     specified frame to another at a specified epoch.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C      FRAMES
C
C$ Keywords
C
C      FRAMES
C
C$ Declarations

      IMPLICIT NONE

      INCLUDE               'zzctr.inc'

      CHARACTER*(*)         FROM
      CHARACTER*(*)         TO
      DOUBLE PRECISION      ET
      DOUBLE PRECISION      ROTATE  ( 3, 3 )

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     FROM       I   Name of the frame to transform from.
C     TO         I   Name of the frame to transform to.
C     ET         I   Epoch of the rotation matrix.
C     ROTATE     O   A rotation matrix
C
C$ Detailed_Input
C
C     FROM        is the name of some reference frame in which
C                 a position vector is known.
C
C     TO          is the name of a reference frame in which it
C                 is desired to represent a position vector.
C
C     ET          is the epoch in ephemeris seconds past the epoch
C                 of J2000 (TDB) at which the position transformation
C                 matrix ROTATE should be evaluated.
C
C$ Detailed_Output
C
C     ROTATE      is the matrix that transforms position vectors from
C                 the reference frame FROM to the frame TO at epoch ET.
C                 If (x, y, z) is a position relative to the frame FROM
C                 then the vector ( x', y', z') is the same position
C                 relative to the frame TO at epoch ET.  Here the
C                 vector ( x', y', z' ) is defined by the equation:
C
C                   -   -       -        -     -  -
C                  | x'  |     |          |   | x  |
C                  | y'  |     |  ROTATE  |   | y  |
C                  | z'  |  =  |          |   | z  |
C                   -   -       -        -     -  -
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1) If sufficient information has not been supplied via loaded
C        SPICE kernels to compute the transformation between the
C        two frames, the error will be diagnosed by a routine
C        in the call tree to this routine.
C
C     2) If either frame FROM or TO is not recognized the error
C        'SPICE(UNKNOWNFRAME)' will be signaled.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     This routine provides the user level interface to computing
C     position transformations from one reference frame to another.
C
C     Note that the reference frames may be inertial or non-inertial.
C     However, the user must take care that sufficient SPICE kernel
C     information is loaded to provide a complete position
C     transformation path from the FROM frame to the TO frame.
C
C$ Examples
C
C     Suppose that you have geodetic coordinates of a station on the
C     surface of the earth and that you need the inertial (J2000)
C     position of this station.  The following code fragment
C     illustrates how to transform the position of the station to a
C     J2000 position.
C
C        CALL BODVRD ( 'EARTH', RADII, 3, N, ABC  )
C
C        EQUATR   =  ABC(1)
C        POLAR    =  ABC(3)
C        F        = (EQUATR - POLAR) / EQUATR
C
C        CALL GEOREC ( LONG, LAT, 0.0D0,  EQUATR, F, EPOS )
C
C        CALL PXFORM ( 'IAU_EARTH', 'J2000', ET,  ROTATE  )
C        CALL MXV    (  ROTATE,      EPOS,   JPOS         )
C
C     The state JPOS is the desired J2000 position of the station.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C      None.
C
C$ Author_and_Institution
C
C     C.H. Acton      (JPL)
C     N.J. Bachman    (JPL)
C     B.V. Semenov    (JPL)
C     W.L. Taber      (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.1.0, 23-SEP-2013 (BVS)
C
C        Updated to save the input frame names and POOL state counters
C        and to do frame name-ID conversions only if the counters have
C        changed.
C
C-    SPICELIB Version 1.0.3, 27-FEB-2008 (BVS)
C
C        Added FRAMES to the Required_Reading section.
C
C-    SPICELIB Version 1.0.2, 23-OCT-2005 (NJB)
C
C        Header example had invalid flattening factor computation;
C        this was corrected.  Reference to BODVAR in header was
C        replaced with reference to BODVRD.
C
C-    SPICELIB Version 1.0.1, 29-JUL-2003 (NJB) (CHA)
C
C        Various header corrections were made.
C
C-    SPICELIB Version 1.0.0, 05-APR-1999 (WLT)
C
C
C-&

C$ Index_Entries
C
C     Find a position transformation matrix
C
C-&

C
C     Spicelib Functions
C
      LOGICAL               RETURN

C
C     Local parameters.
C

C
C     Saved frame name length.
C
      INTEGER               FRNMLN
      PARAMETER           ( FRNMLN = 32 )

C
C     Local Variables.
C
      INTEGER               FCODE
      INTEGER               TCODE

C
C     Saved frame name/ID item declarations.
C
      INTEGER               SVCTR1 ( CTRSIZ )
      CHARACTER*(FRNMLN)    SVFROM
      INTEGER               SVFCOD

      INTEGER               SVCTR2 ( CTRSIZ )
      CHARACTER*(FRNMLN)    SVTO
      INTEGER               SVTCDE

      LOGICAL               FIRST

C
C     Saved frame name/ID items.
C
      SAVE                  SVCTR1
      SAVE                  SVFROM
      SAVE                  SVFCOD

      SAVE                  SVCTR2
      SAVE                  SVTO
      SAVE                  SVTCDE

      SAVE                  FIRST

C
C     Initial values.
C
      DATA                  FIRST   / .TRUE. /

C
C     Standard SPICE error handling.
C
      IF ( RETURN() ) THEN
         RETURN
      END IF

      CALL CHKIN ( 'PXFORM' )

C
C     Initialization.
C
      IF ( FIRST ) THEN

C
C        Initialize counters.
C
         CALL ZZCTRUIN( SVCTR1 )
         CALL ZZCTRUIN( SVCTR2 )

         FIRST = .FALSE.

      END IF

      CALL ZZNAMFRM ( SVCTR1, SVFROM, SVFCOD, FROM, FCODE )
      CALL ZZNAMFRM ( SVCTR2, SVTO,   SVTCDE, TO,   TCODE )

C
C     Only non-zero id-codes are legitimate frame id-codes.  Zero
C     indicates that the frame wasn't recognized.
C
      IF      ( FCODE .NE. 0 .AND. TCODE .NE. 0 ) THEN

         CALL REFCHG ( FCODE, TCODE, ET, ROTATE )

      ELSE IF ( FCODE .EQ. 0 .AND. TCODE .EQ. 0 ) THEN

         CALL SETMSG ( 'Neither of the frames # or # was '
     .   //            'recognized as a known reference frame. ' )
         CALL ERRCH  ( '#', FROM )
         CALL ERRCH  ( '#', TO   )
         CALL SIGERR ( 'SPICE(UNKNOWNFRAME)'  )

      ELSE IF ( FCODE .EQ. 0 ) THEN

         CALL SETMSG ( 'The frame # was not recognized as a '
     .   //            'known reference frame. ' )
         CALL ERRCH  ( '#', FROM )
         CALL SIGERR ( 'SPICE(UNKNOWNFRAME)'  )

      ELSE IF ( TCODE .EQ. 0 ) THEN

         CALL SETMSG ( 'The frame # was not recognized as a '
     .   //            'known reference frame. ' )
         CALL ERRCH  ( '#', TO   )
         CALL SIGERR ( 'SPICE(UNKNOWNFRAME)'  )

      END IF

      CALL CHKOUT ( 'PXFORM' )
      RETURN
      END

ext: f

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
!
! Create_MW_SRF
!
! Program to create netcdf SRF files for the Microwave using data
! from the MW_Sensordata_Define module
!
!
! CREATION HISTORY:
!       Written by:     David Groff, SAIC 01-Jan-2009
!                       david.groff@noaa.gov
! ------------------------------------------------------------------

PROGRAM Create_MW_SRF

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE Type_Kinds
  USE File_Utility
  USE Message_Handler
  USE CRTM_Parameters
  USE MW_SensorData_Define,      ONLY: MW_SensorData_type,      &
                                       Allocate_MW_SensorData,  &
                                       Load_MW_SensorData,      &
                                       Destroy_MW_SensorData
  USE SRF_Define,                ONLY: SRF_type,                &
                                       Integrate_SRF,           &
                                       Allocate_SRF,            &
                                       Destroy_SRF
  USE SRF_netCDF_IO,             ONLY: Create_SRF_netCDF,       &
                                       Write_SRF_netCDF
  USE SensorInfo_IO,             ONLY: Read_SensorInfo
  USE SensorInfo_LinkedList,     ONLY: SensorInfo_List_type,    &
                                       Count_SensorInfo_Nodes,  &
                                       GetFrom_SensorInfo_List
  USE SensorInfo_Define,         ONLY: SensorInfo_type
  USE SensorInfo_Parameters,     ONLY: MICROWAVE_SENSOR

  ! Parameters
  ! ----------
  CHARACTER( * ), PARAMETER :: PROGRAM_RCS_ID = &
  CHARACTER( * ), PARAMETER :: PROGRAM_NAME = 'Create_MW_SRF'

  ! Variables
  ! ---------
  TYPE( MW_SensorData_type ) :: MW_SensorData
  TYPE( SensorInfo_List_type ) :: SensorInfo_List
  TYPE( SRF_type ) :: SRF
  TYPE( SensorInfo_type ) :: SensorInfo
  INTEGER, DIMENSION(:), ALLOCATABLE :: n_SBPoints
  CHARACTER(256)  :: SensorInfo_Filename
  CHARACTER(256)  :: NC_Filename
  CHARACTER(5000) :: Title, Comment
  INTEGER :: n_Frequencies
  INTEGER :: n_Sensors
  INTEGER :: n_Bands
  INTEGER :: n, l, ln
  INTEGER :: Error_Status
  INTEGER :: Version
  REAL(fp) :: df

  ! Program header
  ! --------------
  CALL Program_Message( PROGRAM_NAME,                                            &
                        'Program to create netcdf SRF files for the microwave.', &
                        '$Revision$'                                      )

  ! Get user specified inputs
  ! -------------------------

  ! Set a user specified SensorInfo filename
  WRITE( *,FMT    ='(/5x,"Enter a SensorInfo filename: ")', &
           ADVANCE='NO' )
  READ( *,FMT='(a)' ) SensorInfo_Filename
  SensorInfo_Filename = ADJUSTL(SensorInfo_Filename)

  ! Read SensorInfo file
  Error_Status = Read_SensorInfo( SensorInfo_Filename,  & ! Input
                                  SensorInfo_List       ) ! Output
  IF ( Error_Status /= SUCCESS ) THEN
     CALL Display_Message( PROGRAM_NAME,                     &
                           'Error reading SensorInfo_file',  &
                           Error_Status                      )
    STOP
  END IF

  ! Get number of sensors
  n_Sensors = Count_SensorInfo_Nodes( SensorInfo_List )

  ! Set a user specified version for the SRF files to be created
  WRITE( *,FMT='(/5x,"Default SRF file version is: ",i0, &
               &".  Enter value: ")', &
           ADVANCE='NO' ) SRF%VERSION
  READ( *,* ) Version
  IF ( Version < SRF%VERSION ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Invalid version number specified. Using default.', &
                          INFORMATION )
    Version = SRF%VERSION
  END IF

  ! Get the number of frequency points
  WRITE( *,FMT    ='(/5x,"Enter number of frequencies [256] :")', &
           ADVANCE='NO' )
  READ( *,* ) n_Frequencies

  ! Loop over sensors
  !--------------------------------
  Sensor_Loop: DO n = 1, n_Sensors

    ! Get SensorInfo for sensor n
    Error_Status = GetFrom_SensorInfo_List( SensorInfo_List,     &  ! Input
                                            n,                   &  ! Input
                                            SensorInfo           )  ! Output
    IF ( Error_Status /= SUCCESS ) THEN
       CALL Display_Message( PROGRAM_NAME,                &
                             'Error getting sensor node', &
                             Error_Status                 )
    END IF

    ! Cycle if its not a microwave sensor or if the sensor name is SSMIS
    IF ( .NOT. (SensorInfo%Sensor_Type==MICROWAVE_SENSOR) ) CYCLE Sensor_Loop !.OR. &
    !      TRIM(SensorInfo%Sensor_Name)=='SSMIS' ) CYCLE Sensor_Loop

    ! Name of file to write to for sensor
    NC_Filename= TRIM(SensorInfo%Sensor_Id) // '.srf.nc'

    ! Set the title
    Title = 'Microwave sensor ' // TRIM(SensorInfo%Sensor_Id) // ' Spectral Response Functions'
    Comment = 'The SRF data for this sensor is assumed to be a box car shape.'     // &
              ' This is a temporary state until real data becomes available. The ' // &
              'files now contain 256 points. Further testing will be performed'    // &
              ' to determine if this is an adequate number of points'

    ! Create srf netcdf files for the MW
    Error_Status = Create_SRF_netCDF( NC_Filename                                    ,  & ! Input
                                      MICROWAVE_SENSOR                               ,  & ! Input (MW)
                                      SensorInfo%Sensor_Channel                      ,  & ! Input Channel list
                                      Version          = Version                     ,  & ! Optional Input
                                      Sensor_Id        = TRIM(SensorInfo%Sensor_Id)  ,  & ! Optional Input
                                      WMO_Satellite_Id = SensorInfo%WMO_Satellite_Id ,  & ! Optional Input
                                      WMO_Sensor_Id    = SensorInfo%WMO_Sensor_Id    ,  & ! Optional Input
                                      Title            = Title                       ,  & ! Optional Input
                                      History          = PROGRAM_RCS_ID              ,  & ! Optional Input
                                      Comment          = Comment                        ) ! Optional Input
    IF ( Error_Status /= SUCCESS ) THEN
       CALL Display_Message( PROGRAM_NAME,                      &
                             'Error creating SRF netcdf file',  &
                             Error_Status                       )
    END IF

    ! Fill the MW_SensorData structure
    Error_Status = Load_MW_SensorData( MW_SensorData,                  &
                                       n_Frequencies=n_Frequencies,    &
                                       Sensor_Id=SensorInfo%Sensor_Id  )
    IF ( Error_Status /= SUCCESS ) THEN
       CALL Display_Message( PROGRAM_NAME,                  &
                             'Error loading MW_SensorData', &
                             Error_Status                   )
    END IF

    ! Loop over channels
    ! -------------------------------------------
    Channel_Loop: DO l = 1, SensorInfo%n_Channels

      ! Allocate Sideband points array
      ALLOCATE( n_SBPoints(MW_SensorData%n_Sidebands(l)),  &
                    STAT = Error_Status                    )
      IF ( Error_Status /= SUCCESS ) THEN
         CALL Display_Message( PROGRAM_NAME,                             &
                               'Error allocating for n_SBPoints array',  &
                               Error_Status                              )
      END IF

      ! Set the number of bands for the channel
      IF(MW_SensorData%IF_Band(1,1,l)==ZERO) THEN
        n_Bands = 1
      ELSE
        n_Bands = MW_SensorData%n_Sidebands(l)*2
      END IF

      ! Allocate for SRF structure
      Error_Status = Allocate_SRF( n_Frequencies,  &  ! Input
                                   SRF          ,  &  ! Output
                                   n_Bands         )  ! Optional Output
      IF ( Error_Status /= SUCCESS ) THEN
         CALL Display_Message( PROGRAM_NAME,                &
                               'Error Allocating for SRF',  &
                               Error_Status                 )
      END IF

      ! Assign components of SRF structure
      ! ----------------------------------
      SRF%Channel          = SensorInfo%Sensor_Channel(l)
      SRF%Frequency        = MW_SensorData%Frequency(:,l)
      SRF%Response         = MW_SensorData%Response(:,l)
      SRF%Version          = Version
      SRF%Sensor_Id        = SensorInfo%Sensor_Id
      SRF%WMO_Satellite_ID = SensorInfo%WMO_Satellite_ID
      SRF%WMO_Sensor_ID    = SensorInfo%WMO_Sensor_ID
      SRF%Sensor_Type      = MICROWAVE_SENSOR

      ! Obtain number of points for sidebands
      n_SBPoints(:)=0
      DO ln = 1, MW_SensorData%n_Sidebands( l )
        df = MW_SensorData%IF_Band(2,ln,l) - MW_SensorData%IF_Band(1,ln,l)
        n_SBPoints(ln) = NINT(df/MW_SensorData%Delta_Frequency(l)) + 1
      END DO

      ! assign the number of points to corresponding band and
      ! assign the bounds for the bands
      IF (n_Bands==1) THEN
        SRF%f1_Band(1)=MW_SensorData%Central_Frequency(l)-MW_SensorData%IF_Band(2,1,l)
        SRF%f2_Band(1)=MW_SensorData%Central_Frequency(l)+MW_SensorData%IF_Band(2,1,l)
        SRF%npts_Band(1)=n_Frequencies
      ELSE IF(n_Bands==2) THEN
        SRF%f1_Band(1)=MW_SensorData%Central_Frequency(l)-MW_SensorData%IF_Band(2,1,l)
        SRF%f2_Band(1)=MW_SensorData%Central_Frequency(l)-MW_SensorData%IF_Band(1,1,l)
        SRF%f1_Band(2)=MW_SensorData%Central_Frequency(l)+MW_SensorData%IF_Band(1,1,l)
        SRF%f2_Band(2)=MW_SensorData%Central_Frequency(l)+MW_SensorData%IF_Band(2,1,l)
        SRF%npts_Band(1)=n_SBPoints(1)
        SRF%npts_Band(2)=n_SBPoints(1)
      ELSE IF(n_Bands==4) THEN
        SRF%f1_Band(1)=MW_SensorData%Central_Frequency(l)-MW_SensorData%IF_Band(2,2,l)
        SRF%f2_Band(1)=MW_SensorData%Central_Frequency(l)-MW_SensorData%IF_Band(1,2,l)
        SRF%f1_Band(2)=MW_SensorData%Central_Frequency(l)-MW_SensorData%IF_Band(2,1,l)
        SRF%f2_Band(2)=MW_SensorData%Central_Frequency(l)-MW_SensorData%IF_Band(1,1,l)
        SRF%f1_Band(3)=MW_SensorData%Central_Frequency(l)+MW_SensorData%IF_Band(1,1,l)
        SRF%f2_Band(3)=MW_SensorData%Central_Frequency(l)+MW_SensorData%IF_Band(2,1,l)
        SRF%f1_Band(4)=MW_SensorData%Central_Frequency(l)+MW_SensorData%IF_Band(1,2,l)
        SRF%f2_Band(4)=MW_SensorData%Central_Frequency(l)+MW_SensorData%IF_Band(2,2,l)
        SRF%npts_Band(1)=n_SBPoints(2)
        SRF%npts_Band(2)=n_SBPoints(1)
        SRF%npts_Band(3)=n_SBPoints(1)
        SRF%npts_Band(4)=n_SBPoints(2)
      END IF

      ! Fill the integrated and summation fields
      Error_Status = Integrate_SRF( SRF ) ! In/Output
      IF ( Error_Status /= SUCCESS ) THEN
         CALL Display_Message( PROGRAM_NAME,                                      &
                               'Error Calculating Integrated and Summation SRF',  &
                               Error_Status                                       )
      END IF

      ! Write the SRF data
      Error_Status = Write_SRF_netCDF(  NC_Filename,   &  ! In/Output
                                        SRF            )  ! Input
      IF ( Error_Status /= SUCCESS ) THEN
         CALL Display_Message( PROGRAM_NAME,                                      &
                               'Error writing the SRF data to file',              &
                               Error_Status                                       )
      END IF

      ! Destroy the SRF structure
      Error_Status = Destroy_SRF( SRF )
      IF ( Error_Status /= SUCCESS ) THEN
         CALL Display_Message( PROGRAM_NAME,                                      &
                               'Error destroying SRF structure',                  &
                               Error_Status                                       )
      END IF

      ! Deallocate Sideband points array
      DEALLOCATE( n_SBPoints,          &
                  STAT = Error_Status  )
      IF ( Error_Status /= SUCCESS ) THEN
         CALL Display_Message( PROGRAM_NAME,                               &
                               'Error deallocating for n_SBPoints array',  &
                               Error_Status                                )
      END IF

    END DO Channel_Loop

    ! Destroy MW_SensorData structure
    Error_Status = Destroy_MW_SensorData( MW_SensorData )
    IF ( Error_Status /= SUCCESS ) THEN
         CALL Display_Message( PROGRAM_NAME,                                      &
                               'Error destroying MW_SensorData structure',        &
                               Error_Status                                       )
    END IF

  END DO Sensor_Loop

END PROGRAM Create_MW_SRF

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
! This test code is from test2007_94.f90.
   integer :: a(10)=(/(j,j=1,10)/)
end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
C written by amazing people from PyMC team: https://github.com/pymc-devs/pymc
C Chris Fonnesbeck, Anand Patil, David Huard, John Salvatier
C copied and compiled here, because PyMC is no longer maintained, and PyMC3
C lacks mvhyperg

      DOUBLE PRECISION FUNCTION gammln(xx)
C Returns the value ln[gamma(xx)] for xx > 0.

      DOUBLE PRECISION xx
      INTEGER j
      DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)

C Internal arithmetic will be done in double precision,
C a nicety that you can omit if five-figure accuracy is good enough.

      SAVE cof,stp
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
     +24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
     +-.5395239384953d-5,2.5066282746310005d0/
      x=xx
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*dlog(tmp)-tmp
      ser=1.000000000190015d0
      do j=1,6
         y=y+1.d0
         ser=ser+cof(j)/y
      enddo
      gammln=tmp+dlog(stp*ser/x)
      return
      END

      DOUBLE PRECISION FUNCTION factln(n)
C USES gammln Returns ln(n!).

      INTEGER n
      DOUBLE PRECISION a(100),gammln, pass_val
      DOUBLE PRECISION infinity
      PARAMETER (infinity = 1.7976931348623157d308)

      SAVE a
C Initialize the table to negative values.
      DATA a/100*-1./
      pass_val = n + 1
      if (n.lt.0) then
c        write (*,*) 'negative factorial in factln'
        factln=-infinity
        return
      endif
C In range of the table.
      if (n.le.99) then
C If not already in the table, put it in.
        if (a(n+1).lt.0.) a(n+1)=gammln(pass_val)
        factln=a(n+1)
      else
C Out of range of the table.
        factln=gammln(pass_val)
      endif
      return
      END

      SUBROUTINE mvhyperg(x,color,k,like)

c Multivariate hypergeometric log-likelihood function
c Using the analogy of an urn filled with balls of different colors,
c the mv hypergeometric distribution describes the probability of
c drawing x(i) balls of a given color.
c
c x : (array) Number of draws for each color.
c color : (array) Number of balls of each color.

c Total number of draws = sum(x)
c Total number of balls in the urn = sum(color)

cf2py integer dimension(k),intent(in) :: x,color
cf2py integer intent(hide),depend(x) :: k=len(x)
cf2py double precision intent(out) :: like
cf2py threadsafe

      INTEGER x(k),color(k)
      INTEGER d,total,i,k
      DOUBLE PRECISION like
      DOUBLE PRECISION factln
      DOUBLE PRECISION infinity
      PARAMETER (infinity = 1.7976931348623157d308)

      total = 0
      d = 0
      like = 0.0
      do i=1,k
c Combinations of x balls of color i
        like = like + factln(color(i))-factln(x(i))
     +-factln(color(i)-x(i))
        if ((color(i) .LT. 0.0) .OR. (x(i) .LT. 0.0)) then
          like = -infinity
          RETURN
        endif
        d = d + x(i)
        total = total + color(i)
      enddo
      if (total .LE. 0.0) then
        like = -infinity
        RETURN
      endif
c Combinations of d draws from total
      like = like - (factln(total)-factln(d)-factln(total-d))
      return
      END

max_stars_repo_name: kkedzierska/sonics

max_issues_repo_name: kzkedzierska/sonics

max_forks_repo_name: kzkedzierska/sonics


text:
program SimintTest
use SimintFortran
use iso_c_binding
implicit none
type(c_simint_shell), target :: s_shell(4),p_shell(1)
type(c_simint_multi_shellpair), target :: s2p_msh, s1p_msh

double precision :: alpha(1), coef(1)
double precision :: integrals(1000)
double precision, allocatable :: work(:)

double precision, pointer :: p1(:), p2(:)
integer :: i, ncomputed
integer :: worksize

call simint_init()

worksize = simint_eri_worksize(0, 3)
allocate(work(worksize))

call simint_initialize_shell(s_shell(1))
call simint_initialize_shell(s_shell(2))
call simint_initialize_shell(s_shell(3))
call simint_initialize_shell(s_shell(4))
call simint_initialize_shell(p_shell(1))
call simint_initialize_multi_shellpair(s2p_msh)
call simint_initialize_multi_shellpair(s1p_msh)
! oxygen s
alpha(1)=1d0
coef(1) =1d0
call simint_create_shell(1, 0, 0.0d0, 0.0d0, -0.02d0, &
                          alpha(1), coef(1), s_shell(1))
! oxygen sp
alpha(1)=0.4d0
coef(1) =1d0
call simint_create_shell(1, 0, 0.0d0, 0.0d0, -0.02d0, &
                         alpha(1), coef(1), s_shell(2))
call simint_create_shell(1, 1, 0.0d0, 0.0d0, -0.02d0, &
                        alpha(1), coef(1), p_shell(1))
!hydrogen s
alpha(1)=0.5d0
coef(1) =1d0
call simint_create_shell(1, 0, -0.74d0, 0.0d0, -0.76d0, &
                         alpha(1), coef(1), s_shell(3))
call simint_create_shell(1, 0,  0.74d0, 0.0d0, -0.76d0, &
                         alpha(1), coef(1), s_shell(4))

call C_F_POINTER(s_shell(1)%alpha, p1, shape=[s_shell(1)%nprim])
call C_F_POINTER(s_shell(1)%coef, p2, shape=[s_shell(1)%nprim])
write(*,*) "s_shell1  info"
do i = 1, s_shell(1)%nprim
  write(*,*) p1(i), p2(i)
end do
call C_F_POINTER(s_shell(2)%alpha, p1, shape=[s_shell(2)%nprim])
call C_F_POINTER(s_shell(2)%coef, p2, shape=[s_shell(2)%nprim])
write(*,*) "s_shell1  info"
do i = 1, s_shell(2)%nprim
  write(*,*) p1(i), p2(i)
end do
call C_F_POINTER(p_shell(1)%alpha, p1, shape=[p_shell(1)%nprim])
call C_F_POINTER(p_shell(1)%coef, p2, shape=[p_shell(1)%nprim])
write(*,*) "p_shell  info"
do i = 1, s_shell(1)%nprim
  write(*,*) p1(i), p2(i)
end do

call simint_normalize_shells(4, s_shell)
call simint_normalize_shells(1, p_shell(1))

call simint_create_multi_shellpair(1, s_shell(4), 1, p_shell(1), s2p_msh, 0)
call simint_create_multi_shellpair(1, s_shell(3), 1, p_shell(1), s1p_msh, 0)

call C_F_POINTER(s1p_msh%alpha, p1, shape=[s1p_msh%nprim])
call C_F_POINTER(s1p_msh%prefac, p2, shape=[s1p_msh%nprim])
write(*,*) "s2p Shell Pair info"
do i = 1, s1p_msh%nprim
  write(*,*) p1(i), p2(i)
end do

call C_F_POINTER(s2p_msh%alpha, p1, shape=[s2p_msh%nprim])
call C_F_POINTER(s2p_msh%prefac, p2, shape=[s2p_msh%nprim])
write(*,*) "s2p Shell Pair info"
do i = 1, s2p_msh%nprim
  write(*,*) p1(i), p2(i)
end do

ncomputed = simint_compute_eri(s2p_msh, s1p_msh, 0.0d0, work, integrals)
write(6,*) ' ncomputed ',ncomputed
do i=1,9
   write(*,'(i4,f20.6)') i,integrals(i)
enddo

deallocate(work)
call simint_free_shell(s_shell(1))
call simint_free_shell(s_shell(2))
call simint_free_shell(s_shell(3))
call simint_free_shell(s_shell(4))
call simint_free_shell(p_shell(1))
call simint_free_multi_shellpair(s2p_msh)
call simint_free_multi_shellpair(s1p_msh)

call simint_finalize()

end program

max_stars_repo_name: colleeneb/simint-generator

max_issues_repo_name: colleeneb/simint-generator

max_forks_repo_name: colleeneb/simint-generator


text:
      SUBROUTINE BINEV(IPAIR)
*
*
*       Binary evolution.
*       -----------------
*
      INCLUDE 'common6.h'
      COMMON/BINARY/  CM(4,MMAX),XREL(3,MMAX),VREL(3,MMAX),
     &                HM(MMAX),UM(4,MMAX),UMDOT(4,MMAX),TMDIS(MMAX),
     &                NAMEM(MMAX),NAMEG(MMAX),KSTARM(MMAX),IFLAGM(MMAX)
      REAL*8  M1,M2
      LOGICAL  FIRST
      SAVE  FIRST,IPREV,ISKIP
      DATA  FIRST,IPREV,ISKIP  /.TRUE.,0,0/
*
*
*       Open unit #17 the first time.
      IF (rank.eq.0.and.FIRST) THEN
          OPEN (UNIT=17,STATUS='UNKNOWN',FORM='FORMATTED',FILE='BINEV')
          FIRST = .FALSE.
          WRITE (17,5)
    5     FORMAT ('   TPHYS  NAM1  NAM2  K1  K2  KC   M1   M2     R1',
     &            '    R2   RI    ECC    SEMI   PERIOD  C',/)
      END IF
*
*       Define KS & c.m. indices and set component masses & radii.
      I1 = 2*IPAIR - 1
      I2 = I1 + 1
      I = N + IPAIR
      M1 = BODY(I1)*SMU
      M2 = BODY(I2)*SMU
      R1 = RADIUS(I1)*SU
      R2 = RADIUS(I2)*SU
*
*       Skip identical or alternating sequences (i.e. KSTAR(I) = -1 & 0).
      ISUM = KSTAR(I1) + KSTAR(I2) + KSTAR(I)
      IF (ISUM.EQ.IPREV.OR.ISUM.EQ.ISKIP) GO TO 20
      ISKIP = IPREV
      IPREV = ISUM
*
*       Form basic two-body elements (distinguish mergers & ghosts).
      IF (NAME(I).GT.0.AND.BODY(I).GT.0.0D0) THEN
          SEMI = -0.5*BODY(I)/H(IPAIR)
          BODYC = BODY(I)
          ECC2 = (1.0 - R(IPAIR)/SEMI)**2 + TDOT2(IPAIR)**2/(BODYC*SEMI)
      ELSE
          CALL FINDJ(I1,J,IM)
          IF (J.LE.0) RETURN
          BODYC = CM(1,IM) + CM(2,IM)
          SEMI = -0.5*BODYC/HM(IM)
          RJ = UM(1,IM)**2 + UM(2,IM)**2 + UM(3,IM)**2 + UM(4,IM)**2
          ECC2 = (1.0 - RJ/SEMI)**2
      END IF
*
      ECC = SQRT(ECC2)
      P = DAYS*SEMI*SQRT(ABS(SEMI)/BODYC)
      P = MIN(P,99999.9D0)
      P = MAX(P,-1.0D0)
      A0 = MIN(SEMI*SU,9999.9D0)
      A0 = MAX(A0,-1.0D0)
      RI2 = (X(1,I) - RDENS(1))**2 + (X(2,I) - RDENS(2))**2 +
     &                               (X(3,I) - RDENS(3))**2
      RI = MIN(SQRT(RI2),99.9D0)
      R1 = MIN(R1,999.9D0)
      R2 = MIN(R2,999.9D0)
*
*       Print one line for every new stage (type change of #I1, #I2 or c.m.).
      if(rank.eq.0)then
      WRITE (17,10)  TPHYS, NAME(I1), NAME(I2), KSTAR(I1), KSTAR(I2),
     &               KSTAR(I), M1, M2, R1, R2, RI, ECC, A0, P, IQCOLL
   10 FORMAT (F8.1,2I6,3I4,2F5.1,F7.1,F6.1,F5.1,F7.3,F8.1,F9.1,I3)
      CALL FLUSH(17)
      end if
*
   20 RETURN
*
      END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
      SUBROUTINE SWITCH(Y)
*
*
*       Switching of chain.
*       -------------------
*
      INCLUDE 'commonc.h'
      INCLUDE 'common2.h'
      LOGICAL  KSLOW,KCOLL,Itest
      REAL*8  Y(NMX8),XCNEW(NMX3),KSCH,ksnew(nmx)
      COMMON/SLOW1/   TK2(0:NMX),EJUMP,KSCH(NMX),KSLOW,KCOLL
      INTEGER  IOLD(NMX)
*
*
*       Copy Y-array to COMMON.
      CALL YSAVE(Y)
*
*	First transform to chain coordinates.
      DO I=1,N-1
      L1=3*(I-1)+1
      KS1=4*(I-1)+1
      CALL KSPHYS(Q(KS1),P(KS1),XC(L1),WC(L1))
      END DO
*
      L2=3*(INAME(1)-1)
      DO K=1,3
      X(L2+K)=0.0
      END DO
*
*       Set X for determining new chain indices.
      DO I=1,N-1
      L=3*(I-1)
      L1=L2
      L2=3*(INAME(I+1)-1)
      DO K=1,3
      X(L2+K)=X(L1+K)+XC(L+K)
      END DO
      END DO
*
*	Save the old chain indices.
      DO I=1,N
      IOLD(I)=INAME(I)
      END DO
*
*       Select new indices.
      CALL SELECT
*
*       EXit if chain vectors are unchanged.
      ISW=0
      DO I=1,N
      IF(INAME(I).NE.IOLD(I))ISW=ISW+1
      END DO
      IF(ISW.EQ.0)RETURN
*
*	Transform chain momenta.
      L1=3*(IOLD(1)-1)
      LN=3*(IOLD(N)-1)
      L=3*(N-2)
      DO K=1,3
      PI(L1+K)=-WC(K)
      PI(LN+K)=WC(L+K)
      END DO
      DO I=2,N-1
      L=3*(I-1)
      LI=3*(IOLD(I)-1)
      DO K=1,3
      PI(LI+K)=WC(L+K-3)-WC(L+K)
      END DO
      END DO
      L1=3*(INAME(1)-1)
      LN=3*(INAME(N)-1)
      L=3*(N-2)
      DO K=1,3
      WC(K)=-PI(L1+K)
      WC(L+K)=PI(LN+K)
      END DO
      DO I=2,N-2
      L=3*(I-1)
      LI=3*(INAME(I)-1)
      DO K=1,3
      WC(L+K)=WC(L+K-3)-PI(LI+K)
      END DO
      END DO
*
*       Construct new chain coordinates.
      DO I=1,3*(N-1)
      XCNEW(I)=0.0
      END DO
*       Transformation matrix (old to new) has only coefficients -1, 0 or +1.
      DO ICNEW=1,N-1
*       Find K0 & K1 for IOLD(K0) = INAME(ICNEW) & IOLD(K1) = INAME(ICNEW+1).
      LNEW=3*(ICNEW-1)
      DO I=1,N
      IF(IOLD(I).EQ.INAME(ICNEW))K0=I
      IF(IOLD(I).EQ.INAME(ICNEW+1))K1=I
      END DO
      DO ICOLD=1,N-1
      LOLD=3*(ICOLD-1)
      IF((K1.GT.ICOLD).AND.(K0.LE.ICOLD))THEN
*       Add.
      DO K=1,3
      XCNEW(LNEW+K)=XCNEW(LNEW+K)+XC(LOLD+K)
      END DO
      ELSE IF((K1.LE.ICOLD).AND.(K0.GT.ICOLD))THEN
*	Subtract.
      DO K=1,3
      XCNEW(LNEW+K)=XCNEW(LNEW+K)-XC(LOLD+K)
      END DO
      END IF
      END DO
      END DO
*
*	Perform KS-transformations and update chain coordinates & RINV.
      DO I=1,N-1
      L1=3*(I-1)+1
      KS1=4*(I-1)+1
      CALL PHYSKS(XCNEW(L1),WC(L1),Q(KS1),P(KS1))
      DO K=1,3
      XC(L1+K-1)=XCNEW(L1+K-1)
      END DO
      RINV(I)=1.0/(Q(KS1)**2+Q(KS1+1)**2+Q(KS1+2)**2+Q(KS1+3)**2)
      END DO
*
*       Swop Ksch acording to Iold and Iname.
      do i=1,n-1
          ksnew(i)=1.0
      end do
      do i=1,n-1
          if(ksch(i).ne.1.0d0)then
              i1=Iold(i)
              i2=Iold(i+1)
              do j=1,n-1
              Itest=((i1.eq.iname(j)).and.(i2.eq.iname(j+1))).or.
     &              ((i2.eq.iname(j)).and.(i1.eq.iname(j+1)))
              if(Itest)then
                  ksnew(j)=ksch(i)
              end if
              end do
          end if
      end do
      do i=1,n-1
          ksch(i)=ksnew(i)
      end do
*       Define auxiliary quantities.
      MASS=0.0
      DO I=1,N
      L=3*(I-1)
      MC(I)=M(INAME(I))
      MASS=MASS+MC(I)
      END DO
*
      do i=1,n
          tk1(i)=-1./MC(I)
      end do
      do i=1,n-1
          if(ksch(i).ne.1.0d0)then
              tk1(i)=tk1(i)/Ksch(i)
              tk1(i+1)=tk1(i+1)/Ksch(i)
          end if
      end do
      DO I=1,N-1
          TKK(I)=.5D0*(-tk1(i)-tk1(i+1))
          MKK(I)=MC(I)*MC(I+1)/Ksch(i)
          DO J=I+1,N
              MIJ(I,J)=MC(I)*MC(J)
              MIJ(J,I)=MIJ(I,J)
           END DO
      END DO
      do i=1,n-1
          m12=mc(i)+mc(i+1)
          dt12=0.5d0*(1.d0-1.d0/Ksch(i))/m12
          if(i.gt.1) TKK(i-1)=tkk(i-1)+dt12
          if(i.lt.n-1) TKK(i+1)=tkk(i+1)+dt12
          if(i.gt.1.and.i.lt.n-1) TK2(i)=-2.0d0*dt12
      end do
*       Note: TK2(0) & TK2(N) should be zero but never used.
      TK2(0) = 0.0
      TK2(N) = 0.0
*
*       Copy Y-array from COMMON.
      CALL YCOPY(Y)
*
      RETURN
      END