text:
!***********************************************************************
!                                                                      *
      SUBROUTINE GETCID (isofile, rwffile, idblk)
!                                                                      *
!   Interactively determines the data governing the CI problem.        *
!   iccut is replaced by an array iccutblk(1:nblock)
!                                                                      *
!   Call(s) to: [LIB92]: GETYN, NUCPOT, RADGRD, SETQIC.                *
!               [RCI92]: SETISO, SETRWF.                               *
!                                                                      *
!   Written by Farid A. Parpia            Last revision: 15 Dec 1992   *
!   Block version by Xinghong He          Last revision: 15 Jun 1998   *
!                                                                      *
!***********************************************************************
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:04:58   1/ 3/07
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY: DOUBLE
      USE parameter_def,   ONLY: NNNP
      USE decide_C
      USE def_C,           ONLY: cvac, c, z, accy, nelec, emn
      USE default_C
      USE grid_C
      USE hblock_C
      USE iccu_C
      USE iounit_C
      USE npar_C
      USE npot_C,          ONLY: zz, nnuc
      USE nsmdat_C
      USE orb_C
      USE wave_C
      USE wfac_C
      USE blim_C
      USE where_C
      USE qedcut_C
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      IMPLICIT NONE
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      CHARACTER(LEN=*) :: isofile, rwffile
      CHARACTER(LEN=8), DIMENSION(*) ::  idblk
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
       LOGICAL :: getyn, yes
       INTEGER :: i, jblock, ntmp, np10, j
!-----------------------------------------------------------------------
!
! Open, check, load data from, and close the  .iso  file
! Data loaded are: EMN,Z,/NPAR/,/NSMDAT/ where /.../ means whole
! common block.
!
!      PRINT *, 'Calling SETISO ...'
      CALL SETISO (isofile)
!
! The speed of light, if non-default then spread from node-0
! Quantities to be obtained: C
!
      IF (NDEF .NE. 0) THEN
         WRITE (istde,*) 'Revise the physical speed of light (',CVAC,   &
                           ' in a.u.) ?'
         YES = GETYN ()
         IF (YES) THEN
            WRITE (istde,*) 'Enter the revised value:'
            READ *,C
         ELSE
            C = CVAC
         ENDIF
      ELSE
         C = CVAC
      ENDIF
!
! Treat some CSF's as perturbation ? Broadcasting is used only in
! non-default mode, as the above case for speed of light.
! Quantities to be obtained: LFORDR
!
      IF (NDEF .NE. 0) THEN
         WRITE (istde,*) 'Treat contributions of some CSFs',            &
                           ' as first-order perturbations?'
         LFORDR = GETYN ()
      ELSE
         LFORDR = .FALSE.
      ENDIF

! Get iccutblk() from the user-input

      IF (.NOT. LFORDR) THEN
         !...Default first
         DO i = 1, nblock
           iccutblk(i) = ncfblk(i)
         ENDDO
      ELSE

      ! Let master do the i/o, then broadcast
            WRITE (istde,*) 'There are ', nblock, 'blocks. They are:'
            WRITE (istde,*) '  block     J Parity     No of CSFs'
            DO i = 1, nblock
               WRITE (istde,*) i, idblk(i)(1:5), ncfblk(i)
            ENDDO

            WRITE (istde,*)
            WRITE (istde,*) 'Enter iccut for each block'
            DO jblock = 1, nblock
               WRITE (istde,*) 'Block ', jblock, '   ncf = ',           &
                              ncfblk(jblock)                            &
                      , ' id = ', idblk(jblock)(1:5)
  123          READ (istdi,*) ntmp
               IF (ntmp .GE. 0 .AND. ntmp .LE. ncfblk(jblock)) THEN
                  iccutblk(jblock) = ntmp
               ELSE
                  WRITE (istde,*) 'ICCUT out of range, re-enter:'
                  GOTO 123
               ENDIF
               write(734,*) ntmp,'! ICCUT for block',jblock
            ENDDO
      ENDIF

!*****************************************************************
!
! Pre-run ?
!
!     IF (IPRERUN .EQ. 0) THEN

!        WRITE (istde,*) ' Prerun with limited interaction?'
!        YES = GETYN ()

!        IF (YES) THEN
!           IPRERUN = 1
!           LTRANS = .FALSE.
!           LVP = .FALSE.
!           LNMS = .FALSE.
!           LSMS = .FALSE.
!           LSE = .FALSE.

!           WRITE (istde,*)  ' Give CSL cut'
!           READ *, NCSFPRE
!           WRITE (istde,*)  ' Give coefficient cut for H_0'
!           READ *, COEFFCUT1
!           WRITE (istde,*) ' Give coefficient cut for the transvers'
!    &,                  ' interaction'
!           GOTO 99
!        ENDIF
!     ENDIF
!*****************************************************************
!
! Include transverse ?
! Quantities to be obtained: LTRANS, WFACT
!
      WRITE (istde,*) 'Include contribution of H (Transverse)?'
      LTRANS = GETYN ()
      WRITE (istde,*) 'Modify all transverse photon frequencies?'
      YES = GETYN ()
      IF (YES) THEN
         WRITE (istde,*) 'Enter the scale factor:'
         READ *, WFACT
      ELSE
         WFACT = 1.0D00
      ENDIF
!
! Other interactions ? One logical for each case. Done altogether
!
      WRITE (istde,*) 'Include H (Vacuum Polarisation)?'
      LVP = GETYN ()

      WRITE (istde,*) 'Include H (Normal Mass Shift)?'
      LNMS = GETYN ()

      WRITE (istde,*) 'Include H (Specific Mass Shift)?'
      LSMS = GETYN ()

      WRITE (istde,*) 'Estimate self-energy?'
      LSE = GETYN ()
      IF (LSE.EQV..TRUE.) THEN
         NQEDCUT = 1
         WRITE (istde,*)                                             &
      'Largest n quantum number for including self-energy for orbital'
         WRITE (istde,*) 'n should be less or equal 8'
         READ *, NQEDMAX
      ELSE
         NQEDCUT = 0
      END IF

      IF (LTRANS) THEN
        WRITE(734,'(a)') 'y            ! Contribution of H Transverse?'
      ELSE
        WRITE(734,'(a)') 'n            ! Contribution of H Transverse?'
      END IF
      IF (YES) THEN
        WRITE(734,'(a)') 'y            ! Modify photon frequencies?'
        WRITE(734,*) WFACT,'! Scale factor'
      ELSE
        WRITE(734,'(a)') 'n            ! Modify photon frequencies?'
      END IF
      IF (LVP) THEN
        WRITE(734,'(a)') 'y            ! Vacuum polarization?'
      ELSE
        WRITE(734,'(a)') 'n            ! Vacuum polarization?'
      END IF
      IF (LNMS) THEN
        WRITE(734,'(a)') 'y            ! Normal mass shift?'
      ELSE
        WRITE(734,'(a)') 'n            ! Normal mass shift?'
      END IF
      IF (LSMS) THEN
        WRITE(734,'(a)') 'y            ! Specific mass shift?'
      ELSE
        WRITE(734,'(a)') 'n            ! Specific mass shift?'
      END IF
      IF (LSE) THEN
        WRITE(734,'(a)') 'y            ! Self energy?'
        WRITE(734,*) NQEDMAX, '! Max n for including self energy'
      ELSE
        WRITE(734,'(a)') 'n            ! Self energy?'
      END IF
!
! Parameters controlling the radial grid
!
!cff-cjb grid
!cff use default grid
!
! 99  IF (NPARM .EQ. 0) THEN
!        RNT = EXP (-65.D0/16.D0) / Z
!        H = 0.5D0**4
!        N = MIN (220,NNNP)
!     ELSE
!CFF     .. should be Z-dependent
         RNT = 2.0D-06/Z
         H = 5.D-2
         N = NNNP
!     ENDIF
      HP = 0.D0
!     IF ( NDEF.NE.0) THEN
!        WRITE (istde,*) 'The default radial grid parameters',         &
!                       ' for this case are:'
!        WRITE (istde,*) ' RNT = ',RNT,';'
!        WRITE (istde,*) ' H = ',H,';'
!        WRITE (istde,*) ' HP = ',HP,';'
!        WRITE (istde,*) ' N = ',N,';'
!        WRITE (istde,*) ' revise these values?'
!        YES = GETYN ()
!
!        IF (YES) THEN
!           WRITE (istde,*) 'Enter RNT:'
!           READ *, RNT
!           WRITE (istde,*) 'Enter H:'
!           READ *, H
!           WRITE (istde,*) 'Enter HP:'
!           READ *, HP
!           WRITE (istde,*) 'Enter N:'
!           READ *, N
!        ENDIF
!     ENDIF
!
!   ACCY is an estimate of the accuracy of the numerical procedures
!
      ACCY = H**6
!
!   Set up the coefficients for the numerical procedures
!
      CALL SETQIC
!
!   Generate the radial grid and all associated arrays
!
      CALL RADGRD
!
!   Generate $- r \times V_ (r)$
!
      CALL NUCPOT
!
!   Load the radial wavefunctions
!
      CALL SETRWFA (rwffile)
!
!   Write the basic parameters of the model electron cloud to the
!   .res  file; this is the second record on the file --- the
!   first is the header (see SUBROUTINE SETRES)
!
      WRITE (imcdf) NELEC, NCF, NW, nblock ! ncf is ncftot, from setcsll
!
!   Write the nuclear parameters and $- r \times V_ (r)$
!   to the  .res  file; these are the third, fourth, and fifth
!   records on the file
!
      WRITE (imcdf) Z, EMN
      WRITE (imcdf) NPARM,(PARM(I), I = 1, NPARM)
      WRITE (imcdf) N, (ZZ(I), I = 1, N), NNUC
!
!   Write the physical effects specification to the  .res  file.
!   iccutblk() is now an array of length nblock.
!   This is the sixth record on the file
!
      WRITE (imcdf) C, LFORDR, (ICCUTblk(i), i = 1, nblock),            &
                    LTRANS, WFACT, LVP, LNMS, LSMS
!
!   Write the grid data to the  .res  file; this is the seventh
!   record on the file
!
      NP10 = N + 10
      WRITE (imcdf) RNT, H, HP, (R(I), I = 1, NP10),                    &
                 (RP(I), I = 1, NP10), (RPOR(I), I = 1, NP10) ! (imcdf)
!
!   Write out the interpolated radial wavefunctions; there are
!   2*NW such records; thus the total number of records written
!   at this point is 7+2*NW
!
      DO J = 1, NW
         WRITE (imcdf) E(J), GAMA(J), PZ(J), MF(J)
         WRITE (imcdf) (PF(I,J), I = 1, MF(J)), (QF(I,J), I = 1, MF(J))
      ENDDO
!
      RETURN
      END

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
      subroutine bragg(natms,kmax,cell,xxx,yyy,zzz)

c*********************************************************************
c
c     dl_poly utility routine to calculate the bragg planes in a
c     CONFIG file.
c
c     copyright daresbury laboratory
c     author - w.smith oct 2009
c
c*********************************************************************

      implicit none

      integer i,k,l,m,n,ll,mm,nn,kk,fail,nmin,mmin,natms,kmax,kmax3
      real(8) twopi,peak,stx,sty,stz,det,ckcs,ckss
      real(8) cell(9),xxx(natms),yyy(natms),zzz(natms),rcell(9)

      real(8), allocatable :: ckc(:),cks(:),clm(:),slm(:)
      real(8), allocatable :: emc(:,:),enc(:,:),els(:,:)
      real(8), allocatable :: ems(:,:),ens(:,:),elc(:,:)
      real(8), allocatable :: strk(:)

      data twopi/6.2831853072d0/

      kmax3=((2*kmax+1)**3)/2
      allocate(ckc(natms),cks(natms),clm(natms),slm(natms),stat=fail)
      allocate(emc(natms,0:kmax),enc(natms,0:kmax),stat=fail)
      allocate(els(natms,0:kmax),ems(natms,0:kmax),stat=fail)
      allocate(ens(natms,0:kmax),elc(natms,0:kmax),stat=fail)
      allocate(strk(kmax3),stat=fail)
      if(fail.gt.0)then
        write(*,*)'error - failed array allocation in bragg'
        call exit(0)
      endif

      write(*,"(1x,'bragg plane projection',//,
     x  1x,' l m n   str.fac')")

      peak=1.d0
      open(8,file="BRAGG")

c     invert cell matrix

      call invert(cell,rcell,det)

c     calculate and store exponential factors

      do i=1,natms

        elc(i,0)=1.0d0
        emc(i,0)=1.0d0
        enc(i,0)=1.0d0
        els(i,0)=0.0d0
        ems(i,0)=0.0d0
        ens(i,0)=0.0d0

        stx=rcell(1)*xxx(i)+rcell(4)*yyy(i)+rcell(7)*zzz(i)
        sty=rcell(2)*xxx(i)+rcell(5)*yyy(i)+rcell(8)*zzz(i)
        stz=rcell(3)*xxx(i)+rcell(6)*yyy(i)+rcell(9)*zzz(i)

        elc(i,1)=cos(twopi*stx)
        emc(i,1)=cos(twopi*sty)
        enc(i,1)=cos(twopi*stz)
        els(i,1)=sin(twopi*stx)
        ems(i,1)=sin(twopi*sty)
        ens(i,1)=sin(twopi*stz)

      enddo

      do k=2,kmax

        do i=1,natms

          elc(i,k)=elc(i,k-1)*elc(i,1)-els(i,k-1)*els(i,1)
          emc(i,k)=emc(i,k-1)*emc(i,1)-ems(i,k-1)*ems(i,1)
          enc(i,k)=enc(i,k-1)*enc(i,1)-ens(i,k-1)*ens(i,1)
          els(i,k)=els(i,k-1)*elc(i,1)+elc(i,k-1)*els(i,1)
          ems(i,k)=ems(i,k-1)*emc(i,1)+emc(i,k-1)*ems(i,1)
          ens(i,k)=ens(i,k-1)*enc(i,1)+enc(i,k-1)*ens(i,1)

        enddo

      enddo

c     loop over all k vectors  k=2pi(ll/cl,mm/cl,nn/cl)

      kk=0
      mmin=0
      nmin=1

      do l=0,kmax
        ll=l
        do mm=mmin,kmax
          m=iabs(mm)
          if(mm.ge.0)then
            do i=1,natms
              clm(i)=elc(i,l)*emc(i,m)-els(i,l)*ems(i,m)
              slm(i)=els(i,l)*emc(i,m)+ems(i,m)*elc(i,l)
            enddo
          else
            do i=1,natms
              clm(i)=elc(i,l)*emc(i,m)+els(i,l)*ems(i,m)
              slm(i)=els(i,l)*emc(i,m)-ems(i,m)*elc(i,l)
            enddo
          endif

          do nn=nmin,kmax
            n=iabs(nn)

            if(nn.ge.0)then
              do i=1,natms
                ckc(i)=clm(i)*enc(i,n)-slm(i)*ens(i,n)
                cks(i)=slm(i)*enc(i,n)+clm(i)*ens(i,n)
              enddo
            else
              do i=1,natms
                ckc(i)=clm(i)*enc(i,n)+slm(i)*ens(i,n)
                cks(i)=slm(i)*enc(i,n)-clm(i)*ens(i,n)
              enddo
            endif

            ckcs=0.0d0
            ckss=0.0d0
            do i=1,natms
              ckcs=ckcs+ckc(i)
              ckss=ckss+cks(i)
            enddo
            kk=kk+1
            strk(kk)=(ckcs*ckcs+ckss*ckss)/dble(natms)
            if(strk(kk).gt.peak)
     x        write(*,'(1x,3i4,f10.4)')ll,mm,nn,strk(kk)
            write(8,'(1x,3i4,f10.4)')ll,mm,nn,strk(kk)

          enddo

          nmin=-kmax

        enddo

        mmin=-kmax

      enddo

      deallocate(ckc,cks,clm,slm,strk,stat=fail)
      deallocate(emc,enc,els,ems,ens,elc,stat=fail)
      close(8)

      return
      end
      subroutine invert(a,b,d)

c***********************************************************************
c
c     dl_poly subroutine to invert a 3 * 3 matrix using cofactors
c
c     copyright - daresbury laboratory 1992
c     author    - w. smith       april 1992
c
c***********************************************************************

      implicit none

      real(8) a,b,d,r

      dimension a(9),b(9)

c     calculate adjoint matrix
      b(1)=a(5)*a(9)-a(6)*a(8)
      b(2)=a(3)*a(8)-a(2)*a(9)
      b(3)=a(2)*a(6)-a(3)*a(5)
      b(4)=a(6)*a(7)-a(4)*a(9)
      b(5)=a(1)*a(9)-a(3)*a(7)
      b(6)=a(3)*a(4)-a(1)*a(6)
      b(7)=a(4)*a(8)-a(5)*a(7)
      b(8)=a(2)*a(7)-a(1)*a(8)
      b(9)=a(1)*a(5)-a(2)*a(4)

c     calculate determinant
      d=a(1)*b(1)+a(4)*b(2)+a(7)*b(3)
      r=0.d0
      if(abs(d).gt.0.d0)r=1.d0/d

c     complete inverse matrix
      b(1)=r*b(1)
      b(2)=r*b(2)
      b(3)=r*b(3)
      b(4)=r*b(4)
      b(5)=r*b(5)
      b(6)=r*b(6)
      b(7)=r*b(7)
      b(8)=r*b(8)
      b(9)=r*b(9)

      return
      end subroutine invert

ext: f


text:
Describe Users/UVDRocks here.
Welcome to the Wiki!  Please consider closing this account and opening one using your Importance of using your RealName RealName.  Take a moment and read over the Identity page to see exactly what youre telling other users and readers of the wiki by remaining anonymous.  Particularly when you attack someone who appeared to leave an honest review (even if you disagree with their conclusions), it makes you look bad.  Rather than attacking other users, why not tell readers about your experiences with University Village so they can come to their own conclusions?  Responding with what amounts to repeating your username over and over again doesnt really help your case. Users/TomGarberson

20101116 21:23:06 nbsp Beyond that, you sound like a pretty obvious shill for the company. Users/JoePomidor

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
       SUBROUTINE W3FT40(FLN,GN,PLN,FL,WORK,TRIGS)
C$$$   SUBPROGRAM  DOCUMENTATION  BLOCK
C
C SUBPROGRAM: W3FT40         COMPUTES 2.5 X 2.5 S. HEMI. GRID-SCALER
C   AUTHOR: JONES,R.E.       ORG: W323       DATE: 93-07-23
C
C ABSTRACT: COMPUTES 2.5 X 2.5 S. HEMI. GRID OF 145 X 37 POINTS
C   FROM SPECTRAL COEFFICIENTS IN A RHOMBOIDAL 30 RESOLUTION
C   REPRESENTING A SCALER FIELD.
C
C PROGRAM HISTORY LOG:
C   93-07-23  R.E.JONES   NEW VERSION OF W3FT10, TAKES OUT W3FA12
C                         MAKES PLN 3 DIMENSIONS, PLN IS COMPUTED
C                         ONE TIME IN MAIN PROGRAM, TRADES MEMORY
C                         FOR MORE SPEED. W3FA12 USED 70% OF CPU TIME.
C
C USAGE:  CALL W3FT40(FLN,GN,PLN,FL,WORK,TRIGS)
C
C   INPUT VARIABLES:
C     NAMES  INTERFACE DESCRIPTION OF VARIABLES AND TYPES
C     ------ --------- -----------------------------------------------
C     FLN    ARG LIST  961 COMPLEX COEFF.
C     PLN    ARG LIST  (32,31,37) REAL SPACE WITH LEGENDRE POLYNOMIALS
C                      COMPUTED BY W3FA12.
C     FL     ARG LIST  31 COMPLEX SPACE FOR FOURIER COEFF.
C     WORK   ARG LIST  144 REAL WORK SPACE FOR SUBR. W3FT12
C     TRIGS  ARG LIST  216 PRECOMPUTED TRIG FUNCS. USED
C                      IN W3FT12, COMPUTED BY W3FA13
C
C   OUTPUT VARIABLES:
C     NAMES  INTERFACE DESCRIPTION OF VARIABLES AND TYPES
C     ------ --------- -----------------------------------------------
C     GN     ARG LIST  (145,37) GRID VALUES.
C                      5365 POINT GRID IS TYPE 30 OR 1E O.N. 84
C
C   SUBPROGRAMS CALLED:
C     NAMES                                                   LIBRARY
C     ------------------------------------------------------- --------
C     AIMAG  CMPLX  REAL                                      SYSTEM
C     W3FT12                                                  W3LIB
C
C WARNING: W3FT10 WAS OPTIMIZED TO RUN IN A SMALL AMOUNT OF
C   MEMORY, IT WAS NOT OPTIMIZED FOR SPEED, 70 PERCENT OF THE TIME WAS
C   USED BY SUBROUTINE W3FA12 COMPUTING THE LEGENDRE POLYNOMIALS. SINCE
C   THE LEGENDRE POLYNOMIALS ARE CONSTANT THEY NEED TO BE COMPUTED
C   ONLY ONCE IN A PROGRAM. BY MOVING W3FA12 TO THE MAIN PROGRAM AND
C   COMPUTING PLN AS A (32,31,37) ARRAY AND CHANGING THIS SUBROUTINE
C   TO USE PLN AS A THREE DIMENSION ARRAY THE RUNNING TIME WAS CUT
C   70 PERCENT. ADD FOLLOWING CODE TO MAIN PROGRAM TO COMPUTE EPS, PLN,
C   TRIGS, AND RCOS ONE TIME IN PROGRAM.
C
C       DOUBLE PRECISION EPS(992)    [CHANGE TO REAL ON CRAY]
C       DOUBLE PRECISION COLRA       [CHANGE TO REAL ON CRAY]
C
C       REAL             PLN( 32, 31, 37 )
C       REAL             RCOS(37)
C       REAL             TRIGS(216)
C
C       DATA  PI /3.14159265/
C
C       DRAD = 2.5 * PI / 180.0
C       CALL W3FA11(EPS,30)
C       CALL W3FA13(TRIGS,RCOS)
C       DO LAT = 1,37
C         COLRA = (LAT - 1) * DRAD
C         CALL W3FA12 (PLN(1,1,LAT), COLRA, 30, EPS)
C       END DOC
C
C ATTRIBUTES:
C   LANGUAGE: CRAY CFT77 FORTRAN
C   MACHINE:  CRAY Y-MP8/864, CRAY Y-MP EL2/128
C
C$$$
C
       COMPLEX          FL( 31 )
       COMPLEX          FLN( 31 , 31 )
C
       REAL             GN(145,37)
       REAL             PLN( 32, 31, 37 )
       REAL             TRIGS(216)
       REAL             WORK(144)
C
       SAVE
C
         DO 400 LAT = 1,37
C
           DO 100 L = 1, 31
             FL(L) = (0.,0.)
 100       CONTINUE
C
           DO 300 L = 1, 31
             I = 1
             FL(L) = FL(L)+CMPLX(PLN(I,L,LAT) * REAL(FLN(I,L)) ,
     &       PLN(I,L,LAT) * AIMAG(FLN(I,L)) )
C
             DO  200 I = 2, 30 ,2
               FL(L) = FL(L)-CMPLX(PLN(I,L,LAT) * REAL(FLN(I,L)) ,
     &         PLN(I,L,LAT) * AIMAG(FLN(I,L)) )
               FL(L) = FL(L)+CMPLX(PLN(I+1,L,LAT) * REAL(FLN(I+1,L)),
     &         PLN(I+1,L,LAT) * AIMAG(FLN(I+1,L)))
 200         CONTINUE
C
 300       CONTINUE
C
         CALL W3FT12(FL,WORK,GN(1,LAT ),TRIGS)
 400     CONTINUE
C
         RETURN
       END

max_stars_repo_name: GMAO-SI-Team/jedi-stack

max_issues_repo_name: GMAO-SI-Team/jedi-stack

max_forks_repo_name: GMAO-SI-Team/jedi-stack


text:
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
      SUBROUTINE NSPOPL ( LOGNAM, VERSN )

C
C$ Version
C
C-     Command Loop Configured Version 2.0.0, 10-SEP-1998 (WLT)
C
C         The routine now logs the version of SPICELIB that the
C         program was linked against.
C
C-     Command Loop Configured Version 1.0.0, 3-MAY-1994 (WLT)
C
C         This is the configured version of the Command Loop
C         software as of May 4, 1994
C

      IMPLICIT NONE
      CHARACTER*(*)         LOGNAM
      CHARACTER*(*)         VERSN
C
C     This routine opens the log file that will be used for loging
C     commands.  It should only be called once.  If a log file
C     cannot be opened, the routine will issue a warning message
C     to the default output device.
C
      LOGICAL               HAVE
      INTEGER               POS
      INTEGER               RTRIM
      LOGICAL               EQSTR

      INTEGER               LNGSIZ
      PARAMETER           ( LNGSIZ = 800 )

      INTEGER               LNSIZE
      PARAMETER           ( LNSIZE  = 80 )

      INTEGER               ROOM
      PARAMETER           ( ROOM   = 2 )

      INTEGER               WDSIZE
      PARAMETER           ( WDSIZE = 32 )

      INCLUDE              'commdpar.inc'

      CHARACTER*(FILEN)     LOGFIL
      CHARACTER*(LNGSIZ)    MYERR ( 2 )
      CHARACTER*(LNGSIZ)    REST
      CHARACTER*(LNSIZE)    STYLE
      CHARACTER*(LNSIZE)    ENV
      CHARACTER*(LNSIZE)    TKV
      CHARACTER*(LNSIZE)    ERR

      CHARACTER*(WDSIZE)    ATTR ( ROOM )
      CHARACTER*(WDSIZE)    WAS
      CHARACTER*(WDSIZE)    VALUE
      CHARACTER*(WDSIZE)    IO
      CHARACTER*(WDSIZE)    TIME
      CHARACTER*(WDSIZE)    WARN

      INTEGER               I
      INTEGER               N
      INTEGER               PTR
      INTEGER               START
      LOGICAL               FOUND

      EXTERNAL               NSPWLN
C
C     Empty out the internal error buffers.
C

      MYERR(1) = ' '
      MYERR(2) = ' '

      DO I = 1, ROOM
         ATTR(I) = ' '
      END DO

      CALL NEWFIL ( LOGNAM, 'LOG', LOGFIL )

      IF ( HAVE(MYERR) ) THEN
C
C        See if we can parse the error message as having the
C        string IOSTAT was value imbedded in it.  This isn't
C        pretty, but we can possibly get a better idea of
C        what went wrong this way.
C
         START = POS ( MYERR(1), 'IOSTAT', 1 )

         IF ( START .GT. 0 ) THEN

            REST = MYERR(1)(START:)

            CALL NEXTWD ( REST, IO,    REST )
            CALL NEXTWD ( REST, WAS,   REST )
            CALL NEXTWD ( REST, VALUE, REST )

            IF (      ( EQSTR ( WAS, 'was' ) )
     .          .AND. ( VALUE .NE.   ' '     ) ) THEN

               ERR = ' '
               CALL NPARSI ( VALUE, I, ERR, PTR )

               IF ( ERR .EQ. ' ' ) THEN

                  CALL DCYPHR ( I, FOUND, REST )

                  IF ( FOUND ) THEN
                     MYERR(1)(START:) = REST
                  END IF

               END IF

            END IF
         END IF

         REST = MYERR(1)
         WARN = ' '

         CALL TRNLAT ( 'WARNING',       WARN     )
         CALL TRNLAT ( 'CANNOTOPENLOG', MYERR(2) )

         START = RTRIM ( MYERR(2) )

         CALL PREFIX ( MYERR(2)(1:START), 1, REST )

         STYLE = 'LEFT 1 RIGHT 78 NEWLINE /cr FLAG ' // WARN

         CALL NICEPR_1 ( REST,  STYLE, NSPWLN )
      ELSE
         CALL CURTIM ( TIME          )
         CALL PLTFRM ( ROOM, N, ATTR )
         CALL TKVRSN ( 'TOOLKIT', TKV )
         ENV = ATTR(1)
         CALL SUFFIX ( '---',            1, ENV )
         CALL SUFFIX ( ATTR(2),          1, ENV )
         CALL PREFIX ( 'SPICE Toolkit ', 1, TKV )

         CALL NSPLOG ( ENV,   .TRUE. )
         CALL NSPLOG ( VERSN, .TRUE. )
         CALL NSPLOG ( TKV,   .TRUE. )
         CALL NSPLOG ( TIME,  .TRUE. )
      END IF

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
People say that wiki:WikiPedia:On the Internet, nobody knows youre a dog, but DavisWiki is different because its mostly about things offline in Davis. Offline identities for users here are very common and conventional because many people already know each other offline and it can be unfair to talk about real things and people in Davis, while actually being in Davis, without using our real names.

Besides your offline identity, contributing to the wiki forms an online identity which is a combination of:

  Your Importance of using your RealName username.
  Your User Statistics edits, viewable as your edit history.
  Your comments.
  Any other contributions (positive or negative) that you have made. These may include uploading photographs, adding new relevant pages, etc.

When you choose not to use your offline identity, you make the choice not to be known.  Thats fine; thats your prerogative!  But in so doing, your first statement to everyone here is that you arent willing to engage on a personal level or to stand behind your words.  You create a hurdle that youre going to need to overcome if you want your input to carry the same weight as those of an editor with an established identity.

Unfortunately, the wiki gets its fair share of sockpuppets, Wiki Vandalism vandals, and wiki:wikipedia:Whitewash (censorship) whitewashers.  In order to maintain the integrity of information on the wiki, many editors are on the lookout for such Shadiness Factor shady edits.  And therein lies the catch.  Without any sort of identity your comments carry very little weight and your edits may be interpreted by other editors as vandalism and removed.  The more you participate in the wiki, the stronger your identity will become, for better or for worse.  Since we dont know who you are, it is important to establish a positive identity.  This gives you an enhanced reputation with other Wiki users, and builds trust in what you post.  This is true regardless of your username, but it is particularly important if you choose to remain anonymous.

## Photo public domain from http://commons.wikimedia.org/wiki/File:Nophoto.jpg

Also note that as of January 1, 2011, California has imposed new civil and criminal penalties for impersonating someone online.  It is a criminal act to impersonate another online without their consent for purposes of harming, intimidating, threatening, or defrauding anyonenot just the person being impersonated.  Doing so is punishable by a $1,000 fine and up to a year in jail.  You can also be sued for any harm you cause in the course of the impersonation.  In short, make sure the identity you establish here is your own.


text:
!******************************************************************
! Module used by dcv_deconv_2D.f90, dcv_cgrad_mod.f90 dcv_2D_mod.f90
!
! JLP
! Version 26/03/2001
!******************************************************************
module dcv_2D_data
public
! 16384=128*128
! integer,parameter		:: idim=128*128
 integer,parameter		:: idim=256*256
! integer,parameter		:: idim=512*512
 character(len=5)		:: option
!
! Used by dcv_deconv_2D et dcv_2D_mod:
!
 real,dimension(:),allocatable	:: hh,yy,yy0,yyd,yy_re,yy_im,w0,w1
 real,dimension(:),allocatable	:: hh_re,hh_im,hht_re,hht_im,w_re,w_im
 real				:: alpha
 integer			:: nx,ny,nn
 logical			:: positivity=.false.,simulation=.false.

end module dcv_2D_data

max_stars_repo_name: jlprieur/deconv

max_issues_repo_name: jlprieur/deconv

max_forks_repo_name: jlprieur/deconv


text:
*> \brief \b DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLASV2 + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
*
*       .. Scalar Arguments ..
*       DOUBLE PRECISION   CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLASV2 computes the singular value decomposition of a 2-by-2
*> triangular matrix
*>    [  F   G  ]
*>    [  0   H  ].
*> On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
*> smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
*> right singular vectors for abs(SSMAX), giving the decomposition
*>
*>    [ CSL  SNL ] [  F   G  ] [ CSR -SNR ]  =  [ SSMAX   0   ]
*>    [-SNL  CSL ] [  0   H  ] [ SNR  CSR ]     [  0    SSMIN ].
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] F
*> \verbatim
*>          F is DOUBLE PRECISION
*>          The (1,1) element of the 2-by-2 matrix.
*> \endverbatim
*>
*> \param[in] G
*> \verbatim
*>          G is DOUBLE PRECISION
*>          The (1,2) element of the 2-by-2 matrix.
*> \endverbatim
*>
*> \param[in] H
*> \verbatim
*>          H is DOUBLE PRECISION
*>          The (2,2) element of the 2-by-2 matrix.
*> \endverbatim
*>
*> \param[out] SSMIN
*> \verbatim
*>          SSMIN is DOUBLE PRECISION
*>          abs(SSMIN) is the smaller singular value.
*> \endverbatim
*>
*> \param[out] SSMAX
*> \verbatim
*>          SSMAX is DOUBLE PRECISION
*>          abs(SSMAX) is the larger singular value.
*> \endverbatim
*>
*> \param[out] SNL
*> \verbatim
*>          SNL is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[out] CSL
*> \verbatim
*>          CSL is DOUBLE PRECISION
*>          The vector (CSL, SNL) is a unit left singular vector for the
*>          singular value abs(SSMAX).
*> \endverbatim
*>
*> \param[out] SNR
*> \verbatim
*>          SNR is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[out] CSR
*> \verbatim
*>          CSR is DOUBLE PRECISION
*>          The vector (CSR, SNR) is a unit right singular vector for the
*>          singular value abs(SSMAX).
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date September 2012
*
*> \ingroup auxOTHERauxiliary
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  Any input parameter may be aliased with any output parameter.
*>
*>  Barring over/underflow and assuming a guard digit in subtraction, all
*>  output quantities are correct to within a few units in the last
*>  place (ulps).
*>
*>  In IEEE arithmetic, the code works correctly if one matrix element is
*>  infinite.
*>
*>  Overflow will not occur unless the largest singular value itself
*>  overflows or is within a few ulps of overflow. (On machines with
*>  partial overflow, like the Cray, overflow may occur if the largest
*>  singular value is within a factor of 2 of overflow.)
*>
*>  Underflow is harmless if underflow is gradual. Otherwise, results
*>  may correspond to a matrix modified by perturbations of size near
*>  the underflow threshold.
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
*
*  -- LAPACK auxiliary routine (version 3.4.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     September 2012
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION   CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
*     ..
*
* =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0D0 )
      DOUBLE PRECISION   HALF
      PARAMETER          ( HALF = 0.5D0 )
      DOUBLE PRECISION   ONE
      PARAMETER          ( ONE = 1.0D0 )
      DOUBLE PRECISION   TWO
      PARAMETER          ( TWO = 2.0D0 )
      DOUBLE PRECISION   FOUR
      PARAMETER          ( FOUR = 4.0D0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            GASMAL, SWAP
      INTEGER            PMAX
      DOUBLE PRECISION   A, CLT, CRT, D, FA, FT, GA, GT, HA, HT, L, M,
     $                   MM, R, S, SLT, SRT, T, TEMP, TSIGN, TT
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, SIGN, SQRT
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           DLAMCH
*     ..
*     .. Executable Statements ..
*
      FT = F
      FA = ABS( FT )
      HT = H
      HA = ABS( H )
*
*     PMAX points to the maximum absolute element of matrix
*       PMAX = 1 if F largest in absolute values
*       PMAX = 2 if G largest in absolute values
*       PMAX = 3 if H largest in absolute values
*
      PMAX = 1
      SWAP = ( HA.GT.FA )
      IF( SWAP ) THEN
         PMAX = 3
         TEMP = FT
         FT = HT
         HT = TEMP
         TEMP = FA
         FA = HA
         HA = TEMP
*
*        Now FA .ge. HA
*
      END IF
      GT = G
      GA = ABS( GT )
      IF( GA.EQ.ZERO ) THEN
*
*        Diagonal matrix
*
         SSMIN = HA
         SSMAX = FA
         CLT = ONE
         CRT = ONE
         SLT = ZERO
         SRT = ZERO
      ELSE
         GASMAL = .TRUE.
         IF( GA.GT.FA ) THEN
            PMAX = 2
            IF( ( FA / GA ).LT.DLAMCH( 'EPS' ) ) THEN
*
*              Case of very large GA
*
               GASMAL = .FALSE.
               SSMAX = GA
               IF( HA.GT.ONE ) THEN
                  SSMIN = FA / ( GA / HA )
               ELSE
                  SSMIN = ( FA / GA )*HA
               END IF
               CLT = ONE
               SLT = HT / GT
               SRT = ONE
               CRT = FT / GT
            END IF
         END IF
         IF( GASMAL ) THEN
*
*           Normal case
*
            D = FA - HA
            IF( D.EQ.FA ) THEN
*
*              Copes with infinite F or H
*
               L = ONE
            ELSE
               L = D / FA
            END IF
*
*           Note that 0 .le. L .le. 1
*
            M = GT / FT
*
*           Note that abs(M) .le. 1/macheps
*
            T = TWO - L
*
*           Note that T .ge. 1
*
            MM = M*M
            TT = T*T
            S = SQRT( TT+MM )
*
*           Note that 1 .le. S .le. 1 + 1/macheps
*
            IF( L.EQ.ZERO ) THEN
               R = ABS( M )
            ELSE
               R = SQRT( L*L+MM )
            END IF
*
*           Note that 0 .le. R .le. 1 + 1/macheps
*
            A = HALF*( S+R )
*
*           Note that 1 .le. A .le. 1 + abs(M)
*
            SSMIN = HA / A
            SSMAX = FA*A
            IF( MM.EQ.ZERO ) THEN
*
*              Note that M is very tiny
*
               IF( L.EQ.ZERO ) THEN
                  T = SIGN( TWO, FT )*SIGN( ONE, GT )
               ELSE
                  T = GT / SIGN( D, FT ) + M / T
               END IF
            ELSE
               T = ( M / ( S+T )+M / ( R+L ) )*( ONE+A )
            END IF
            L = SQRT( T*T+FOUR )
            CRT = TWO / L
            SRT = T / L
            CLT = ( CRT+SRT*M ) / A
            SLT = ( HT / FT )*SRT / A
         END IF
      END IF
      IF( SWAP ) THEN
         CSL = SRT
         SNL = CRT
         CSR = SLT
         SNR = CLT
      ELSE
         CSL = CLT
         SNL = SLT
         CSR = CRT
         SNR = SRT
      END IF
*
*     Correct signs of SSMAX and SSMIN
*
      IF( PMAX.EQ.1 )
     $   TSIGN = SIGN( ONE, CSR )*SIGN( ONE, CSL )*SIGN( ONE, F )
      IF( PMAX.EQ.2 )
     $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, CSL )*SIGN( ONE, G )
      IF( PMAX.EQ.3 )
     $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, SNL )*SIGN( ONE, H )
      SSMAX = SIGN( SSMAX, TSIGN )
      SSMIN = SIGN( SSMIN, TSIGN*SIGN( ONE, F )*SIGN( ONE, H ) )
      RETURN
*
*     End of DLASV2
*
      END

max_stars_repo_name: pereirak/mpc

max_issues_repo_name: pereirak/mpc

max_forks_repo_name: pereirak/mpc


text:
!
!     CalculiX - A 3-dimensional finite element program
!     Copyright (C) 1998-2021 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     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, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine compdt(nk,dt,nshcon,shcon,vold,ntmat_,
     &     iponoel,inoel,dtimef,ielmat,dh,cocon,
     &     ncocon,ithermal,mi,
     &     vcon,compressible,tincf,ierr,ifreesurface,dgravity,iit)
!
!     - determine the time step for each node (stored in field dt
!     and the minimum value across all nodes (dtimef)
!
      implicit none
!
      integer nk,i,iponoel(*),inoel(2,*),index,nelem,ithermal(*),
     &     mi(*),compressible,ierr,ifreesurface,iit,
     &     nshcon(*),ntmat_,ielmat(mi(3),*),imat,ncocon(2,*),
     &     iflag
!
      real*8 dtimef,dt(*),dvi,r,cp,rho,shcon(0:3,ntmat_,*),tincf,
     &     vold(0:mi(2),*),temp,vel,dtcon,dtmed,dtimefold,
     &     dh(*),cocon(0:6,ntmat_,*),dtthd,cond,c,
     &     vcon(nk,0:mi(2)),dgravity
!
      data iflag /3/
!
!     determining the time increment dt for each node.
!
!     edge nodes (fields iponoel and inoel are determined in precfd.f)
!
      dtimefold=dtimef
      dtimef=1.d30
!
      do i=1,nk
        index=iponoel(i)
        if(index.le.0) cycle
!
!     look at an element belonging to the edge node
!
        nelem=inoel(1,index)
!
!     determining the time increment
!
        imat=ielmat(1,nelem)
        temp=vold(0,i)
!
!     convective time step (dtcon)
!
        vel=dsqrt(vold(1,i)**2+vold(2,i)**2+vold(3,i)**2)
!
        if(compressible.eq.1) then
          if(ifreesurface.eq.0) then
!
!     gas: speed of sound = dsqrt(kappa*r*temp)
!
            rho=vcon(i,4)
            call materialdata_cp(imat,ntmat_,temp,shcon,nshcon,cp)
            r=shcon(3,1,imat)
            dtcon=dh(i)/(dsqrt(cp*r*temp/(cp-r))+vel)
            c=cp*r*temp/(cp-r)
          else
!
!     shallow water: speed of sound = dsqrt(g*depth)
!
            dtcon=dh(i)/(dsqrt(dgravity*vcon(i,4))+vel)
          endif
        else
!
!     liquid
!
          rho=vcon(i,4)
          if(vel.lt.1.d-10) vel=1.d-10
          dtcon=dh(i)/vel
        endif
        dt(i)=dtcon
!
!     mechanical diffusion time step (dtmed)
!
        call materialdata_dvifem(imat,ntmat_,temp,shcon,nshcon,dvi)
        if(dvi.gt.1.d-20) then
          dtmed=dh(i)*dh(i)*rho/(2.d0*dvi)
          dt(i)=dtcon*dtmed/(dtcon+dtmed)
        endif
!
!     thermal diffusion time step (dtthd)
!
        if(ithermal(1).gt.1) then
          call materialdata_cond(imat,ntmat_,temp,cocon,ncocon,
     &         cond)
          call materialdata_cp(imat,ntmat_,temp,shcon,nshcon,cp)
          if(cond.gt.1.d-20) then
            dtthd=dh(i)*dh(i)*rho*cp/(2.d0*cond)
            dt(i)=(dt(i)*dtthd)/(dt(i)+dtthd)
          endif
        endif
!
!     safety factor (default for compressible fluids: 1.25,
!     default for incompressible fluids: 1.00;
!     can be changed by user: fifth entry underneath *CFD)
!
        dt(i)=dt(i)/tincf
c        dt(i)=dt(i)/10.d0
!
        if(dt(i).lt.dtimef) then
          dtimef=dt(i)
        endif
!
      enddo
!
!     increased damping for incompressible fluids
!     the use of an internal (for the damping) and an external
!     (for the time derivative) time step stems from Zienkiewicz,
!     Taylor and Nithiarasu, The Finite Element Method for Fluid
!     Dynamics, 6th edition, p94 bottom and p 95 top.
!
      write(*,*)
      write(*,*) 'dtimef= ',dtimef
      write(*,*)
      if(dtimef.le.0.d0) then
        write(*,*) '*ERROR in compdt;'
        write(*,*) '       negative time increment;'
        write(*,*) '       the solution diverged'
        ierr=1
      elseif(((iit.gt.1).and.(dtimef.lt.dtimefold/100.d0))) then
        write(*,*) '*ERROR in compdt;'
        write(*,*) '       strongly decreasing time increment;'
        write(*,*) '       the solution diverged'
        ierr=1
      endif
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
c--- These are the functions similar to equations (A.8)-(A.14)
c--- but which describe the case of a pseudoscalar (CP odd) Higgs
c--- These results are adapted from Spira et al., hep-ph/9504378.

      double complex function ehsva4_odd(s,t,u)
C     ehsv:EqnA.8
      implicit none
      double precision s,t,u
      double complex ehsvb4_odd
      ehsva4_odd=ehsvb4_odd(s,t,u)+ehsvb4_odd(u,s,t)+ehsvb4_odd(t,u,s)
      return
      end

      double complex function ehsva2_odd(s,t,u)
C     ehsv:EqnA.9
      implicit none
      double precision s,t,u
      double complex ehsvb2_odd
      ehsva2_odd=ehsvb2_odd(s,t,u)+ehsvb2_odd(s,u,t)
      return
      end

      double complex function ehsvb4_odd(s,t,u)
      implicit none
C     ehsv:EqnA.10
      include 'masses.f'
      double precision hmass2,s,t,u
      double complex w2,w3
      hmass2=s+t+u
      ehsvb4_odd=mbsq/hmass2*(w2(hmass2)-w2(s)-w3(s,t,u,hmass2))/6d0
      return
      end

      double complex function ehsvb2_odd(s,t,u)
C     ehsv:EqnA.11
      implicit none
      include 'masses.f'
      double precision hmass2,s,t,u
      double complex w2,w3
      hmass2=s+t+u
      ehsvb2_odd=mbsq/hmass2**2*(
     .                       -s*w3(s,t,u,hmass2)+s/2d0*w3(t,s,u,hmass2)
     .                       +2d0*s*(0.75d0-u/(s+u))*w2(hmass2)
     .                       -s/2d0*w2(s)
     .                       -s*(1d0-2d0*u/(s+u))*w2(t))/6d0
      return
      end

      double complex function ehsva5_odd(s,t,u)
C     ehsv:EqnA.14
      implicit none
      include 'masses.f'
      double precision hmass2,s,t,u
      double complex w2
      hmass2=s+t+u
      ehsva5_odd=mbsq/hmass2*2d0/3d0*(w2(s)-w2(hmass2))
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
      SUBROUTINE DCINFO (IH, NH, IBGELM, NELMS, ISTWRD, NWORDS, IBBYTE)
C
      PARAMETER ( NCHARS = 5 )
C
      COMMON /WORDS/ NCMW, NCPW, IWDS(8)
C
C     COMMON / DCDBUG / LDEBUG, MUNIT                                   D
C     LOGICAL LDEBUG                                                    D
C
      LOGICAL LRPEAT, LDELTA, LSIGDT
C
      INTEGER IH(*), ISTBYT(NCHARS)
      INTEGER NH, IBGELM, NWORDS, IBBYTE
C
C     INTEGER*3 NELMS                                                   H
      INTEGER*4 NELMS                                                   MLu
C
      CHARACTER*1 CVAR, CHEAD(5)
C
C --- TURN ON THE DEBUG OPTION IF ISTAT IS SET TO -999,
C --- THEN RESET ISTAT TO 0
C
C
      ISTAT = 0
C
C --- CHECK TO SEE IF NUMBER OF INPUT ELEMENTS IS LESS THAN OR EQUAL
C --- TO ZERO. IF SO, PRINT ERROR MESSAGE AND RETURN.
C
      IF ( NELMS .LE. 0 ) THEN
         ISTAT = -1201
C        IF (LDEBUG) THEN                                               D
C           WRITE (MUNIT,*) 'NUMBER OF INPUT ELEMENTS LESS THAN 1 '     D
C           WRITE (MUNIT,*) 'NELMS:    ', NELMS, ' ISTAT: ', ISTAT      D
C        ENDIF                                                          D
         GOTO 900
      ENDIF
C
C --- CHECK TO SEE IF THE HEADER HAS BEEN PASSED. IF NH IS LESS THAN
C --- OR EQUAL TO ZERO, PRINT ERROR MESSAGE AND RETURN.
C
      IF ( NH .LE. 0 ) THEN
         ISTAT = -1202
C        IF (LDEBUG) THEN                                               D
C           WRITE (MUNIT,*) 'NUMBER OF ELEMENTS IN HEADER IS LESS '     D
C           WRITE (MUNIT,*) 'THAN OR EQUAL TO ZERO.'                    D
C           WRITE (MUNIT,*) 'NH: ', NH, ' ISTAT: ', ISTAT               D
C        ENDIF                                                          D
         GOTO 900
      ENDIF
C
C --- IF SPECIFIED BEGINNING ELEMENT IS LESS THAN OR EQUAL TO ZERO,
C --- PRINT ERROR MESSAGE AND RETURN TO CALLING PROGRAM.
C
      IF ( IBGELM .LE. 0 ) THEN
         ISTAT = -1203
C        IF (LDEBUG) THEN                                               D
C           WRITE (MUNIT,*) 'NUMBER OF INPUT ELEMENTS LESS THAN 1 '     D
C           WRITE (MUNIT,*) 'NELMS:    ', NELMS, ' ISTAT: ', ISTAT      D
C        ENDIF                                                          D
         GOTO 900
      ENDIF
C
C --- GET THE SECOND BYTE OF THE COMPRESSION HEADER AND VERIFY
C --- THAT THE COMPRESSION SCHEME VERSION NUMBER IS VALID.
C
      CALL GETHOL ( IH, 2, IVERS )
C
      IF ( IVERS .NE. 16 ) THEN
         ISTAT = -1204
C        IF (LDEBUG) THEN                                               D
C           WRITE (MUNIT,*) 'INVALID DATA COMPRESSION VERSION'          D
C           WRITE (MUNIT,*) 'IVERS:  ', IVERS, ' ISTAT: ', ISTAT        D
C        ENDIF                                                          D
         GOTO 900
      ENDIF
C
C --- INITIALIZE THE CHEAD ARRAY WHICH REPRESENTS ASCII CHARACTERS
C --- FOUND IN THE DATA COMPRESSION HEADER.
C
      CHEAD(1) = 'C'
      CHEAD(2) = 'D'
      CHEAD(3) = 'S'
      CHEAD(4) = 'R'
      CHEAD(5) = 'E'
C
C --- CHECK THE FIRST BYTE OF EACH SECTION OF THE COMPRESSION HEADER
C --- FOR THE CORRECT ASCII CHARACTER. IF WRONG CHARACTER IS FOUND,
C --- PRINT ERROR MESSAGE AND RETURN TO CALLING PROGRAM.
C
      ISTBYT(1) = 1
C
      DO 100 I=1,NCHARS
C
      CALL GETHOL ( IH, ISTBYT(I), NBYTES )
      CVAR(1:1) = CHAR( NBYTES )
C
      IF ( CVAR(1:1) .NE. CHEAD(I) ) THEN
         ISTAT = -1205
C        IF (LDEBUG) THEN                                               D
C           WRITE (MUNIT,*) 'INVALID DATA COMPRESSION HEADER'           D
C           WRITE (MUNIT,*) 'CVAR: ', CVAR, ' ISTAT: ', ISTAT           D
C        ENDIF                                                          D
         GOTO 900
      ENDIF
C
      IBYT = ISTBYT(I) + 2
      CALL GETHOL ( IH, IBYT, NBYTES )
      ISTBYT(I+1) = ISTBYT(I) + NBYTES
C
  100 CONTINUE
C
C --- HEADER IS INTACT, SO DETERMINE WHICH COMPRESSION SCHEMES ARE USED
C
      LRPEAT = .FALSE.
      LDELTA = .FALSE.
      LSIGDT = .FALSE.
C
      CALL GETHOL ( IH, 4, ISKM )
      JSKM = ISKM
C
      DO 200 I=7,0,-1
         IBIT = JSKM/2**I
         IF (IBIT .NE. 0) THEN
            JSKM = MOD (JSKM, 2**I)
            IF ( JSKM .EQ. 0 ) LRPEAT = .TRUE.
            IF ( JSKM .EQ. 1 ) LDELTA = .TRUE.
            IF ( JSKM .EQ. 2 ) LSIGDT = .TRUE.
         ENDIF
  200 CONTINUE
C
C --- IF INFORMATION ON THE REPEAT SCHEME IS REQUESTED, THE ENTIRE
C --- ARRAY MUST BE UNCOMPRESSED BEFORE DETERMINING LOCATION OF IBGELM.
C --- THEREFORE PRINT MESSAGE THAT THIS IS NOT IMPLEMENTED.
C
      IF ( LRPEAT ) THEN
         ISTAT = -1206
C        IF (LDEBUG) THEN                                               D
C           WRITE (MUNIT,*) 'DCINFO IS NOT IMPLEMENTED FOR '            D
C           WRITE (MUNIT,*) 'THE REPEAT COUNT SCHEME.      '            D
C           WRITE (MUNIT,*) 'ISKM: ', ISKM, ' ISTAT: ', ISTAT           D
C        ENDIF                                                          D
         GOTO 900
      ENDIF
C
C --- IF REQUESTING INFORMATION ON AN INVALID COMPRESSION TECHNIQUE,
C --- PRINT ERROR MESSAGE AND RETURN TO THE CALLING PROGRAM.
C
      IF ( .NOT. LDELTA .AND. .NOT. LSIGDT ) THEN
         ISTAT = -1207
C        IF (LDEBUG) THEN                                               D
C           WRITE (MUNIT,*) 'INVALID DATA COMPRESSION METHOD'           D
C           WRITE (MUNIT,*) 'ISKM: ', ISKM, ' ISTAT: ', ISTAT           D
C        ENDIF                                                          D
         GOTO 900
      ENDIF
C
C --- GET ISIZE FOR THE DATA COMPRESSION SCHEME USED
C
      IF (LDELTA) THEN
         CALL GETHOL ( IH, ISTBYT(2)+3, ISIZE )
      ELSEIF (LSIGDT) THEN
         CALL GETHOL ( IH, ISTBYT(3)+3, ISIZE )
      ENDIF
C
C --- USING ISIZE AND IBGELM, CALCULATE THE STARTING INTEGER WORD
C --- WHERE COMPRESSION INFORMATION LIVES.
C
      IBBYTE = IBGELM*ISIZE - ISIZE + 1
C
      IREMAN = MOD(IBBYTE,NCPW)
C
      IF ( IREMAN .EQ. 0 ) THEN
         ISTWRD = IBBYTE/NCPW
      ELSE
         ISTWRD = IBBYTE/NCPW + 1
      ENDIF
C
C --- NOW CALCULATE THE NUMBER OF INTEGER WORDS REQUIRED TO CONTAIN
C --- NELMS OF COMPRESSED DATA FOR THE DELTA OR SIGNIFICANT DIGIT
C --- SCHEMES
C
      NBYTES = NELMS*ISIZE
C
      IENBYT = IBBYTE + NBYTES - 1
C
      IREMAN = MOD(IENBYT,NCPW)
C
      IF ( IREMAN .EQ. 0 ) THEN
         IENWRD = IENBYT/NCPW
      ELSE
         IENWRD = IENBYT/NCPW + 1
      ENDIF
C
      NWORDS = IENWRD - ISTWRD + 1
C
C --- RETURN TO CALLING PROGRAM WITH LOCATION AND QUANTITY OF BYTES USED
C --- TO STORE THE COMPRESSED DATA.
C
  900 RETURN
      END

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


text:
      real*4 function t(h, q, z, lamb)
c...... Returns age of universe at redshift z
      implicit none
      real*4 h, q, z, a, b, c, hh0, cosh, x, lamb
c... for lambda
      real*4 aa, bb, epsr, s0, s, funq, omega0
      integer npts
      logical success
      external funq
      common /cosm/ omega0

c	H = Ho in km/sec/Mpc
c	Q = qo  (if problems with qo = 0, try 0.0001)

      a(q,z) = (sqrt(1. + ((2. * q) * z)) / (1. - (2. * q))) / (1. + z)
      b(q) = q / (abs((2. * q) - 1.) ** 1.5)
      c(q,z) = ((1. - (q * (1. - z))) / q) / (1. + z)

      cosh(x) = alog(x + sqrt((x ** 2) - 1.))
      hh0 = h * 0.001022
c     in (billion years)**(-1)

      if (lamb .ne. 0.0) then
         	omega0 = (2. * (q + 1.)) / 3.
         	aa = 0.
         	bb = 1. / (1. + z)
	        success=.false.
		s0=1.e-10
		npts=0
		do while (.not.success)
		  npts=npts+1
		  callmidpnt(funq,aa,bb,s,npts)
		  epsr=abs(s-s0)/s0
		  if (epsr.lt.1.e-4) then
			success=.true.
		  else
			s0=s
		  endif
		enddo
		t=s
      else if (q .eq. 0.0) then
      	t = 1. / (1. + z)
      else if (q .lt. 0.5) then
      	t = a(q,z) - (b(q) * cosh(c(q,z)))
      else if (q .eq. 0.5) then
      	t = (2. / 3.) / ((1. + z) ** 1.5)
      else
      	t = a(q,z) + (b(q) * cos(c(q,z)))
      end if

      t = t / hh0

      return
      end
c=======================================================================
      real*4 function funq(x)
c.... For non-zero cosmological constant
      real*4 x, omega0, omegainv
      common /cosm/ omega0
      omegainv = 1. / omega0
      funq = sqrt(omegainv) / (x*sqrt((omegainv-1.)+(1./(x**3.))))
      return
      end

max_stars_repo_name: jchavesmontero/galaxpy

max_issues_repo_name: jchavesmontero/galaxpy

max_forks_repo_name: jchavesmontero/galaxpy


text:
PROGRAM band_plot

  implicit none

  integer,parameter :: u1=5,u2=110,u3=120,max_loop=100000
  integer,parameter :: jyuu=10
  integer,allocatable :: mbi(:),maxloc10(:,:,:),indx(:)
  integer,allocatable :: itmp(:),itmp1(:),itmp2(:),jyuu_tmp(:,:)
  integer,allocatable :: itmp4(:),jtmp4(:),itmp3(:),jtmp3(:)
  integer,allocatable :: iiii(:),jjjj(:),kkkk(:),itp(:),jtp(:)
  integer :: n,mb_max,mb_ovlp,mb_min,mbv,lll,ll,loop,ll0,lll0
  integer :: i,j,k,l,nbk,mbc,mb1,mb2,mb,jj0,j1,i1,ndata
  integer :: i0,iwork(8),jwork(8),idummy,idummy2,mbk,iblk,nblk
  integer :: jj,jj1,ii1,i_fixed,l0,ierr,m1,m0,m0_max,m1_max
  real(8),parameter :: HT=27.2116d0
  real(8),allocatable :: eval(:,:),esp(:,:,:),dk_bz(:,:),dval(:,:)
  real(8),allocatable :: sqovlp10(:,:,:),rtmp3(:),rtmp4(:)
  real(8) :: rdummy(3),rr0,dummy,evb,ecb,r1,bb(3,3),r

!  open(u1,file="band_eigv_tot",status='old')
!  open(u2,file="band_ovlp_tot",status='old',form="unformatted")
!  open(u3,file="band_dedk_tot",status='old')

  read(u1,*) bb(1:3,1)
  read(u1,*) bb(1:3,2)
  read(u1,*) bb(1:3,3)

  mb_max=0
  do loop=1,max_loop
     read(u1,*) mbv,mb
     write(*,*) "mb,mbv =",mb,mbv
     mb_max=max(mb_max,mb)
     do i=1,mb
        read(u1,*)
     end do
     read(u1,*,END=10)
     backspace(u1)
  end do
  stop "stop!!! too much data"
10 continue
  nbk=loop

  write(*,*) "mb_max =",mb_max
  write(*,*) "nbk    =",nbk

  allocate( mbi(nbk)            ) ; mbi=0
  allocate( dk_bz(3,nbk)        ) ; dk_bz=0.d0
  allocate( eval(mb_max,nbk)    ) ; eval=0.d0
  allocate( dval(mb_max,nbk)    ) ; dval=0.d0
  allocate( esp(nbk,mb_max,nbk) ) ; esp=0.d0

  rewind u1
  read(u1,*)
  read(u1,*)
  read(u1,*)
  do k=1,nbk
     read(u1,*) idummy,mbi(k),dk_bz(1:3,k)
     do i=1,mbi(k)
        read(u1,*) idummy,eval(i,k),dummy
     end do
  end do

  goto 9
  m1_max=0
  m0_max=0
  rewind u2
  do k=2,nbk
     read(u2) m1,m0
     write(*,*) k,m1,m0
     if ( m1 > m1_max .or. m0 > m0_max ) then
        if ( allocated(maxloc10) ) deallocate(maxloc10)
        if ( allocated(sqovlp10) ) deallocate(sqovlp10)
        m1_max=max(m1,m1_max)
        m0_max=max(m0,m0_max)
        allocate( maxloc10(m1_max,m0_max,nbk) )
        allocate( sqovlp10(m1_max,m0_max,nbk) )
     end if
     read(u2) maxloc10(1:m1,1:m0,1)
     read(u2) sqovlp10(1:m1,1:m0,1)
  end do
  maxloc10(:,:,:)=0
  sqovlp10(:,:,:)=0.d0
  rewind u2
  do k=2,nbk
     read(u2) m1,m0
     write(*,*) k,m1,m0
     read(u2) maxloc10(1:m1,1:m0,k)
     read(u2) sqovlp10(1:m1,1:m0,k)
  end do

  rewind u3
  do k=1,nbk
     read(u3,*)
     do i=1,mbi(k)
        read(u3,*) idummy,dummy,rdummy(1:2),dval(i,k)
     end do
  end do
9 continue
!  close(u1)
!  close(u2)
!  close(u3)

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

!  call plot_eval(mb1,mb2,ecb,HT)

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

!  do k=1,1
!     do i=1,mbi(k)
!        write(*,'(1x,2i6,f10.5,5(1x,i4,f10.5))') k,i,eval(i,k)
!     end do
!  end do
!  do k=2,nbk
!     do i=1,mbi(k)
!        write(*,'(1x,2i6,f10.5,5(1x,i4,f10.5))') &
!             k,i,eval(i,k),(maxloc10(i,j,k),sqovlp10(i,j,k),j=1,3)
!     end do
!  end do
!  stop

  rewind 10
  n=minval(mbi)
  do i=1,maxval(mbi)
  do k=1,nbk
     if ( eval(i,k) == 0.d0 ) cycle
     write(10,'(1x,i8,50f10.5)') k,eval(i,k)
  end do
  write(10,*)
  write(10,*)
  end do
  stop

  esp(:,:,:)=1.d10

  iblk=0

  do mbk=nbk,1,-1

     iblk=iblk+1
     write(*,*) "mbk=",mbk

     do i=1,mbi(mbk)
        if ( eval(i,mbk) == 1.d10 ) cycle
        esp(mbk,i,iblk)=eval(i,mbk)
        eval(i,mbk)=1.d10
     end do

     do i=1,mbi(mbk)
        i0=i
        write(*,*) "i0=",mbk,i0
        do k=mbk,2,-1
           j=maxloc10(i0,1,k)
           r=sqovlp10(i0,1,k)
           if ( r < 0.3d0 ) exit
           if ( eval(j,k-1) == 1.d10 ) exit
           write(*,*) "i0=",k-1,j,r
           esp(k-1,i,iblk)=eval(j,k-1)
           eval(j,k-1)=1.d10
           i0=j
        end do
     end do

  end do ! mbk

  nblk=iblk

  j=0
  do iblk=1,nblk
     do i=1,maxval(mbi)
        if ( all(esp(:,i,iblk)==1.d10) ) cycle
        do k=1,nbk
           if ( esp(k,i,iblk) == 1.d10 ) cycle
           write(10,'(1x,i8,50f10.5)') k,esp(k,i,iblk)
        end do
        write(10,*)
        write(10,*)
        j=j+1
        write(10,'("#",i6)') j
     end do
  end do

  write(*,*) count(eval==1.d10),count(eval==0.d0),count(eval/=0.d0),size(eval)

  rewind 11
  write(11,'(1x,"plot \")')
  write(11,'(1x,a16,i4," w lp \")') ' "fort.10" index',0
  do n=1,j
     write(11,'(1x,a16,i4," w lp \")') ',"fort.10" index',n
  end do
  stop

END PROGRAM band_plot

max_stars_repo_name: j-iwata/RSDFT

max_issues_repo_name: j-iwata/RSDFT

max_forks_repo_name: j-iwata/RSDFT


text:
: zcount ( zaddr - addr n )   dup dup if  65535 0 scan drop over - then ;
: zlength ( zaddr - n )   zcount nip ;
: zplace ( from n to - )   tuck over + >r  cmove  0 r> c! ;
: zappend ( from n to - )   zcount + zplace ;
[undefined] third [if] : third  >r over r> swap ; [then]
[undefined] @+ [if] : @+  dup @ swap cell+ swap ; [then]
: u+  rot + swap ;  \ "under plus"
: ?lit  state @ if postpone literal then ;
: do postpone ?do ; immediate
: for  s" 0 ?do" evaluate ; immediate
: /allot  here over allot swap erase ;
: allotment  here swap /allot ;
: move,   here over allot swap move ;
: h?  @ h. ;
: reclaim  h ! ;
: ]#  ] postpone literal ;
: <<  s" lshift" evaluate ; immediate
: >>  s" rshift" evaluate ; immediate
: bit  dup constant  1 lshift ;
: clamp  ( n low high - n ) -rot max min ;
: and!  dup >r @ and r> ! ;
: or!   dup >r @ or r> ! ;
: xor!   dup >r @ xor r> ! ;
: not!   >r invert r> and! ;
: @!  dup @ >r ! r> ;
: bounds  over + swap ;
: lastbody  last @ name> >body ;
: ccount  dup c@ 1 u+ ;
: .name  dup if body> >name ccount type space else . then ;
: $=  compare 0= ;

: count  dup @ cell u+ ;
: string,  dup , move, ;
: place  2dup ! cell+ swap move ;
: append  2dup 2>r count + swap move 2r> +! ;
: count!  ! ;
: count+!  +! ;
: ,"  [char] " parse string, ;
: included  2dup cr ." [Include] " type included ;

\ lo and hi are inclusive
: inrange ( n lo hi - flag )  over - >r - r> #1 + u< ;

: ifill  ( addr count val - )  -rot  0 do  over !+  loop  2drop ;
: ierase   0 ifill ;
: imove  ( from to count - )  cells move ;
: time?  ( xt - ) ucounter 2>r  execute  ucounter 2r> d-  d>s  . ;

: kbytes  #1024 * ;
: megs    #1048576 * ;
: udup  over swap ;
: 2,  swap , , ;
: 3,  rot , swap , , ;
: 4,  2swap swap , , swap , , ;
: :make  :noname  postpone [  [compile] is  ] ;
: reverse   ( ... count - ... ) 1 + 1 ?do i 1 - roll loop ;
: ;then  s" exit then" evaluate ; immediate
: free  dup 0= if ;then free ;

\ Random numbers
0 VALUE seed
: /rnd  ucounter drop to seed ;  /rnd
: random ( - u ) seed $107465 * $234567 + DUP TO seed ;
: rnd ( n - 0..n-1 ) random um* nip ;

\ readability helper: slang words.  callable once then they self-destruct.
: ?compile  state @ if compile, else execute then ;
: does-slang  does> dup @ ?compile  0 swap body> >name c! ;
: :slang  ( - <name> <code> ; )   create immediate here 0 , does-slang  :noname swap ! ;

\ vocabulary helpers
: define
    >in @
    exists if >in ! also ' execute definitions exit then  \ already defined
    dup >in !  vocabulary
    >in !  also ' execute definitions ;
: using  only forth definitions also ;
vocabulary internal

\ on-stack vector stuff (roger)
: 2!  swap over cell+ ! ! ;
: 2@  dup @ swap cell+ @ ;
: 2+!  swap over cell+ +! +! ;
: 3@  dup @ swap cell+ dup @ swap cell+ @ ;
: 4@  dup @ swap cell+ dup @ swap cell+ dup @ swap cell+ @ ;
: 3!  dup >r  2 cells + !  r> 2! ;
: 4!  dup >r  2 cells + 2! r> 2! ;
: 2min  rot min >r min r> ;
: 2max  rot max >r max r> ;
: 2+  rot + >r + r> ;
: 2-  rot swap - >r - r> ;
: 2negate  negate swap negate swap ;
: 2clamp  ( x y lowx lowy highx highy - x y ) 2>r 2max 2r> 2min ;

\ Word tools
: defined ( - c-addr 0 | xt -1 | - xt 1 )  bl word find ;
: exists ( - flag )   defined 0 <> nip ;

\ compile and exec
: :now  :noname  [char] ; parse evaluate  postpone ;  execute ;

include afkit/ans/depend.f

defer alert  ( a c - )
:make alert  type true abort ;

max_stars_repo_name: RogerLevy/RamenEngine

max_issues_repo_name: Alim-Oezdemir/RamenEngine

max_forks_repo_name: Alim-Oezdemir/RamenEngine


text:
!------------ cutting routine for DM
      subroutine dm_cuts(p,failed)
      implicit none
      include 'constants.f'
!      include 'dm_params.f'
      include 'jetlabel.f'
      include 'leptcuts.f'
      double precision p(mxpart,4)
      logical failed,first
      double precision met_min,met_calc
      data first /.true./
      save first
      double precision dphi_jj
      double precision phi_jjmax

!======= default is to pass

      failed=.false.

      met_min=misspt
      phi_jjmax=2.5d0
!      pt_jmin=100d0

!====== missing ET cuts
!====== met p_T(3,4) always

      met_calc=(p(3,1)+p(4,1))**2+(p(3,2)+p(4,2))**2
      met_calc=dsqrt(met_calc)

      if(met_calc.lt.met_min) failed=.true.

      if(first) then
         first=.false.
         write(6,*) '************ DM CUTS ****************'
         write(6,55) 'cutting MET at ',met_min
!         write(6,55) 'cutting jet pt at ',pt_jmin
         write(6,55) 'Delta Phi <  ',phi_jjmax
         write(6,*) '*************************************'
      endif

      if(jets.eq.2) then
!--- cut on phi_jj
        dphi_jj=
     .   (p(5,1)*p(6,1)+p(5,2)*p(6,2))
     .   /dsqrt((p(5,1)**2+p(5,2)**2)
     .         *(p(6,1)**2+p(6,2)**2))
      if (dphi_jj .lt. -0.999999999D0) dphi_jj=-1d0
      dphi_jj=dacos(dphi_jj)
      if(dphi_jj.gt.phi_jjmax) failed=.true.
      endif

 55   format(1x,a20,f8.3)

      return
      end


text:
Fire University(F.U.) is a weekly community run fire event that teaches how to safely fire dance/spin for free.  There are not formal classes, but rather every student is a teacher.  It typically starts around 9pm and goes until 11pm.  On clear nights it is held at E Street Plaza and during rainy nights it moves to Central Park.

Fire University was started in 2006.  It is not a performance troupe and does not accept monetary donations.  Everyone is welcome to come to watch and learn as long as they are willing to follow some basic community principles.  Alcohol containers are not tolerated at F.U.

To spin fire one must dress only in natural fibers and never nylon or polyester.  There is a minimum proficiency that should be achieved with a non fire prop before anyone lights up.  When going out to spin fire, ensure that you have a designated fire safety ready to use a duvetyne to extinguish any unwanted fire.  F.U. is a self managing community and everyone present is asked to confront anyone who is acting in an unsafe manner or blatantly breaking laws.

The fuel of choice at F.U. is Coleman Stove Fuel aka White Gas available at ACE Hardware.  There is a centralized fuel dump with a communal dip can.  This lighted area is where all fuel should be kept and dispensed.

If you wish to learn, there will be people willing to teach you a variety of props.

Dont forget your fuel!!!


text:
      SUBROUTINE taiutc ( IY, IM, ID, FD, DELTAT, J )
*
*  - - - - - - - -
*   taiutc
*  - - - - - - - -
*
*  This routine is a MODIFIED version of the  the International Astronomical Union's
*  SOFA (Standards of Fundamental Astronomy) software collection.
*
*  Status:  user-replaceable support routine.
*
*  Given:
*     IY       i     UTC:  year (Notes 1 and 2)
*     IM       i           month (Note 2)
*     ID       i           day (Notes 2 and 3)
*     FD       d           fraction of day (Note 4)
*
*  Returned:
*     DELTAT   d     TAI minus UTC, seconds
*     J        i     status (Note 5):
*                       1 = dubious year (Note 1)
*                       0 = OK
*                      -1 = bad year
*                      -2 = bad month
*                      -3 = bad day (Note 3)
*                      -4 = bad fraction (Note 4)
*                      -5 = internal error (Note 5)
*
*  Notes:
*
*  1) UTC began at 1960 January 1.0 (JD 2436934.5) and it is improper
*     to call the routine with an earlier date.  If this is attempted,
*     zero is returned together with a warning status.
*
*  1a)In this verion of the routine, dates before the first entry in the
*     leap.seconds file (January 1st 1972) cannot be used.
*
*     Because leap seconds cannot, in principle, be predicted in
*     advance, a reliable check for dates beyond the valid range is
*     impossible.  To guard against gross errors, a year five or more
*     after the release year of the present routine (see parameter IYV)
*     is considered dubious.  In this case a warning status is returned
*     but the result is computed in the normal way.
*
*     For both too-early and too-late years, the warning status is J=+1.
*     This is distinct from the error status J=-1, which signifies a
*     year so early that JD could not be computed.
*
*  2) If the specified date is for a day which ends with a leap second,
*     the TAI-UTC value returned is for the period leading up to the
*     leap second.  If the date is for a day which begins as a leap
*     second ends, the TAI-UTC returned is for the period following the
*     leap second.
*
*  3) The day number must be in the normal calendar range, for example
*     1 through 30 for April.  The "almanac" convention of allowing
*     such dates as January 0 and December 32 is not supported in this
*     routine, in order to avoid confusion near leap seconds.
*
*  4) The fraction of day is used only for dates before the introduction
*     of leap seconds, the first of which occurred at the end of 1971.
*     It is tested for validity (0 to 1 is the valid range) even if not
*     used;  if invalid, zero is used and status J=-4 is returned.  For
*     many applications, setting FD to zero is acceptable;  the
*     resulting error is always less than 3 ms (and occurs only
*     pre-1972).
*
*  5) The status value returned in the case where there are multiple
*     errors refers to the first error detected.  For example, if the
*     month and day are 13 and 32 respectively, J=-2 (bad month) will be
*     returned.  The "internal error" status refers to a case that is
*     impossible but causes some compilers to issue a warning.
*
*  6) In cases where a valid result is not available, zero is returned.
*
*  Called:
*     iau_CAL2JD   Gregorian calendar to JD
*
*  This revision:  2019 July 9th
*
*-----------------------------------------------------------------------

      USE mdl_param
      IMPLICIT NONE

      INTEGER IY, IM, ID
      DOUBLE PRECISION FD, DELTAT
      INTEGER J

*  Miscellaneous local variables
      LOGICAL MORE
      INTEGER JS, M, N, IS
      DOUBLE PRECISION DA, DJM0, DJM

* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

*  Initialize the result to zero and the status to OK.
      DA = 0D0
      JS = 0

*  If invalid fraction of a day, set error status and give up.
      IF ( FD.LT.0D0 .OR. FD.GT.1D0 ) THEN
         JS = -4
         GO TO 9000
      END IF

*  Convert the date into an MJD.
      CALL iau_CAL2JD ( IY, IM, ID, DJM0, DJM, JS )

*  If invalid year, month, or day, give up.
      IF ( JS .LT. 0 ) GO TO 9000

*  If pre-UTC year, set warning status and give up.
      IF ( IY .LT. IDAT(1,1) ) THEN
         JS = 1
         GO TO 9000
      END IF

*  If suspiciously late year, set warning status but proceed.
      IF ( IY .GT. IDAT(1,NDAT)+3 ) JS = 1

*  Combine year and month.
      M = 12*IY+IM

*  Find the most recent table entry.
      IS = 0
      MORE = .TRUE.
      DO 1 N=NDAT,1,-1
         IF ( MORE ) THEN
            IS = N
            MORE = M .LT. ( 12*IDAT(1,N) + IDAT(2,N) )
         END IF
 1    CONTINUE

*  Prevent underflow warnings.
      IF ( IS .LT. 1 ) THEN
         JS = -5
         GO TO 9000
      END IF

*  Get the Delta(AT).
      DA = DATS(IS)

*  Return the Delta(AT) value and the status.
 9000 CONTINUE
      DELTAT = DA
      J = JS
!      PRINT *,'TAI-UTC: ',DELTAT

*  Finished.

      END

ext: for


text:
 Thank you!

http://www.kickstarter.com/projects/455852114/localwikibringcollaborativelocalmediatoevery/pledge/new Pledge today and help us bring collaborative, communityowned media to every community!

Our fundraiser drive is allornothing  if we dont meet our goal, and cant really do this right, you wont have to pay anything!

Weve also got a number of great Wiki Pledge Drive 2010/Pledge Rewards rewards lined up, so if you pledge  in addition to helping us make this project truly successful  youll get a little (or a lot) of something!

   Users/JasonMoore Jason Moore
   vinh pham
   Debra Van Sickle
   Users/AndrewLeonard Andrew Leonard
   Stephen Fairclough
   Users/GrahamFreeman Graham Freeman
   Users/RodneyBlackwell Rodney Blackwell
   Users/HannahMoore Hannah Moore
   Gary Neustrom
   Users/ArlenAbraham Arlen Abraham
   Lizzy
   Michael Andersen
   Users/PaulIvanov Paul Ivanov
   Users/KenjiYamada Kenji Yamada
   Users/AndrewChen Andrew Chen
   Jim Ratliff
   Kimberlie Kranich
   Users/TerryCliff Terry Cliff
   Matt Lewis
   Nicole Aptekar
   Iosh Heller
   Users/WilliamLewis William Lewis
   Users/Babhari Bahareh Abhari
   Users/SarahEdwards Sarah and Users/JabberWokky Evan Edwards
   Yawen Kurniawan
   Rebekka and Users/TomGarberson Tom Garberson
   Erin Digitale
   Users/KristenBirdsall Kristen Birdsall
   Users/JasonAller Jason Aller
   Bliss Bruen
   Lauren Giles
   Users/CarlosBarahona Carlos Barahona
   Users/ScottMeehleib Scott Meehleib
   Users/JoePomidor Joe Pomidor
   Users/NickSchmalenberger Nick Schmalenberger
   Geoff
   Users/ARWENNHOLD Adam Wennhold
   Users/RobRoy Rob Roy
   Users/BrentLaabs Brent Laabs
   Ryan Deussing
   Sue Hjerpe
   Marilyn Macie
   David Vernon
   Sheryl Schaffner
   Lee Yeung
   Users/DagonJones Dagon Jones
   Users/AmyZimmerman Amy Zimmerman
   Users/VictorZapanta Victor Zapanta
   Ben Antonius
   Users/KalenGallagher Kalen Gallagher
   Roberta
   Ryan Massie
   Users/BlairBaldwin Blair Baldwin
   Users/MichellePalmer Michelle Palmer
   Users/CharlesMcLaughlin Charles McLaughlin
   Eric Toscano
   York Ng
   Adam Dewit
   Don and Vi Schaffner
   Users/JeffTolentino Jeff Tolentino
   Matteo Mattone
   Benny Brown
   Users/csroberson Charles Roberson
   Kika
   Luke Grecki
   Stephanie Thrasher
   Dale Kunce
   Users/TheShah Smit Shah
   Eugene Eric Kim
   Nihar Shah
   Clara
   Maia
   Ryan Wold
   Tara Nonymous
   Rick
   Davis Bike Collective
   Users/AmeliaCarlson Amelia Gray
   Emily Goyins
   Kathleen
   Users/ScottLay Scott Lay
   http://www.wikihow.com/ WikiHow
   Jose D Lopez
   Cheyenne Cary
   Laura Wynholds
   Users/Darach Darach Miller
   Users/SallyHensel Sally Hensel
   Aaron Christopher Cohen
   Ink Monkey Graphics, who donated a beautiful 6 x 3 DavisWiki banner for Farmers Market tabling
   Users/BradBenedict Brad Benedict
   Yumiko
   Kurt Vaughn
   Wrye Sententia
   David Dickson
   Dina Biscotti and Siva Gunda
   Nico Linesch
   Craig Wells
   Reragatz
   June Huang
   Olivia Youngers
   Robert Barat
   Tree
   Sharon Shewmake
   Megan Clements
   Mikkel Marius Winther
   Peter Boulay (RocWiki)
   Sharon
   Besha Gray
   Monica Vaz
   Bonnie Bilger
   Users/EmilyTung Emily Tung
   Craig Blomberg
   The Domes
   Kao Saechao
   Christine Hjelmfelt
   Ippie
   Dave Menninger
   Users/KevinC Kevin Chuc
   Ian Leary
   Ray Porfiri
   Aezlo
   Simon
   Overland Music Distributors
   Dan Ryan
   Cobalt Salon & Gallery
   Paul Lopes
   Jesse Woo
   Efrem Rensi
   Brian
   Angela Clark / Street Corner
   Princess N Aaron Curtin / Apex Cycles
   Users/NaBarry Nicholas Barry
   Liz Ernst
   Sarah White
   Sarah Manley
   Ashley Smith
   Yogesh Simmhan
   Doug Beeferman
   Bernadette Balics / Ecological Landscape Design
   Joady
   Peter Holmes à Court
   Users/TedBuehler Ted Buehler
   Divonne Holmes à Court
   Paige Craig
   Davis Food Coop
   Redwood Barn Nursery
   Mal Jago
   Diane Cary
   Tony Finnis
   Erie Vitiello
   Wyatt Pearsal
   Tobin Fricke
   Christopher Stollery
   Dorai Thodla
   Christopher Michael
   Barbara Mckee
   Anthony Bruce
   Hendrik
   Kao Saechao
   Elle Cayabyab Gitlin
   Virginia Matzek
   Kevin Hawins
   Tori
   Alex Tafreshi
   Daniel Tancredi
   Kris Bulman
   Rosanna Yau
   Kathy
   Jasmine Carlisle
   Andrew Frishman
   Anna Kostecke
   Michael Milman
   Mort Clark
   Sharon Zimmerman
   Bernard Kreissman
   Jasmin Horn
   Michael Trice
   Bev Sykes
   Scott Weintraub
   Stuart Thomas
   Nancy Wilson
   Joyce Maguire
   Abbie
   Michael Johnston
   Mike Adams
   Jessica Winter
   Steve
   Craig Brozinsky
   Michelle Kang
   Gordon Windeyer
   Sunny Grace
   Max Yaffa
   William Gray
   Simon
   Nirgunan Tiruchelvam
   Con Karson
   Ben Gray
   Geoff Ball
   Ivan
   Mark Flitcroft
   A Roach
   Anton Pound
   Derrick
   Misty
   Gavin Payne
   Katharine
   Sandy
   Eric Blosch
   Paul Siegel
   Mehdi Vahab
   Zac Morris
   Rachel Rogers
   Users/EmmyMelton Emmy Melton
   Roberta M.
   Diane C Murphey
   Jonathan London
   Nori
   Matt Quan
   Users/itayneta Itay Neta / Davis Community Acupuncture Clinic
   Cressida Downey
   Lois
   Zach Short
   Yell
   Joseph
   Hank Childs
   Jodi Morrison
   Saroja Subrahmanyan
   Tim Kerbavaz
   Harry Spanglet
   Kevin Liu
   Sam Logan
   Andrew Logan
   Monique Logan
   Emily Logan
   Jamie
   Terry Bellerby
   Stu Nash
   Gabriel Goldstein
   Meredith Clark
   Joel
   Users/MichelleAccurso Michelle Accurso
   Francoise
   Jessica Lao
   Erik Magnuson
   Allison
   Julian Hofer
   Eric Lee
   Tom Denton
   Malcolm
   Tomas Z
   Linda
   Amy Catania
   Charlotte Goodenough
   Chris Lambertus
   Rob Gysel
   Derek Downey
   MeiLyn Peterson
   Mike Bellamy
   Craig Brady
   Nathanael Bettridge
   Bronte Brodie

Want to help out with our pledge drive?  Email philipn@gmail.com for more info!!

20100910 16:35:59 nbsp I just wanted to commend everybody whos donated, especially in the last week.  10 days ago, the $25K goal looked like a distant star, but to see us with just less than $1000 to go and three more days to do it in I feel like its actually going to happen!  Yay for Daviswiki! Users/jefftolentino
  I keep thinking next time, lets not do an allornothing drive again, but it sure does feel good as the finish line approaches.  XNot there yet, however!X jw
   Well, the allornothing approach does add some extra motivation.  In this case, it seemed to work though.  I dont think there would have been as much response, by both donors and volunteers, if it there was a chance of losing it.  I think it was a good thing, maybe slightly stressful to some extent, but still a good thing.
    Woohoo!!! cp

20100910 19:35:57 nbsp WOOWOOOO!!! Users/PhilipNeustrom

20100910 21:01:48 nbsp Gratz gang.. will try to get some $$ to you, just got hit with another child support whammy.  ouch. Users/WesP

20100911 09:04:52 nbsp WOOOHOOOOO!!! Users/jefftolentino

20100911 13:59:47 nbsp Image(JWoo.jpg, noborder, thumbnail, 350) Users/JabberWokky
 Best. Picture.  Ever.  cp


module clock_class
    implicit none
    private

    type :: a_clock
        real :: current ! current time
        real :: total   ! total integration time
        real :: step    ! integration step
        real :: timer   ! cyclic time count for printing out
        real :: lap     ! cyclic period for printing out
    contains
        procedure, public :: init
        procedure, public :: next_step
        procedure, public :: timer_has_elapsed
        procedure, public :: time_is_not_over
        procedure, public :: get_step
    end type a_clock

    type(a_clock), public :: clock

contains

    !---------------------------------------------------------------------------
    ! - read-in the integration step, the cyclic lap period, and the total time;
    ! - set the current time to zero;
    ! - set the cyclic time count to the cyclic lap period.
    subroutine init(self)
        class(a_clock), intent(inout) :: self
        read *, self%step
        read *, self%lap
        read *, self%total
        self%current = 0.0
        self%timer   = self%lap
        return
    end subroutine init

    !---------------------------------------------------------------------------
    ! - increment the current time by one step;
    ! - decremet the cyclic timer by one step.
    subroutine next_step(self)
        class(a_clock), intent(inout) :: self
        self%current = self%current + self%step
        self%timer   = self%timer   - self%step
        return
    end subroutine next_step

    !---------------------------------------------------------------------------
    logical function timer_has_elapsed(self) result(has_elapsed)
        class(a_clock), intent(inout) :: self
        has_elapsed = ( self%timer <= 0.0 )
        if (has_elapsed) self%timer = self%timer + self%lap
        return
    end function timer_has_elapsed

    !---------------------------------------------------------------------------
    logical function time_is_not_over(self)
        class(a_clock), intent(in) :: self
        time_is_not_over = ( self%current <= self%total )
        return
    end function time_is_not_over

    !---------------------------------------------------------------------------
    real function get_step(self) result(step)
        class(a_clock), intent(in) :: self
        step = self%step
        return
    end function get_step

end module clock_class