text:
      SUBROUTINE PLATTR (TSTYPE)

C  PLATTR tests polyline attributes, either individually, using the
C  default values for the other attributes, or in combination, using
C  non-default values for the others.

C  Input parameter:
C    TSTYPE : 'i' for individual, 'c' for combination

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

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

      COMMON /DIALOG/ DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT,
     2                DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS
      INTEGER         DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT
      REAL            DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS

C aspect source
C                bundled     individual
      INTEGER    PBUNDL,     PINDIV
      PARAMETER (PBUNDL = 0, PINDIV = 1)

C linetype
      INTEGER    PLSOLI,     PLDASH,     PLDOT,     PLDASD
      PARAMETER (PLSOLI = 1, PLDASH = 2, PLDOT = 3, PLDASD = 4)

C type of returned value
      INTEGER    PSET,     PREALI
      PARAMETER (PSET = 0, PREALI = 1)

C colour model
      INTEGER    PRGB,     PCIE,     PHSV,     PHLS
      PARAMETER (PRGB = 1, PCIE = 2, PHSV = 3, PHLS = 4)

      INTEGER    IMPMAX
      PARAMETER (IMPMAX = 600)

      INTEGER    PICSTR,TXCI, NUMLT,NUMLW, STDDX,REGDX,IMPDX, LDX
      INTEGER    THISLT, LAVSLT(4), LAVRLT(100), LAVILT(IMPMAX)
      INTEGER    ITRIM, MINILT,MAXILT, UNSUDX,LUNSUP(10), NUMDIS
      INTEGER    LINDIS(20), CDIM,CSIZ, SZCOLT, NUMLIN, NUMPAS
      INTEGER    LNCOL(20), IX, VISDX,VISCOL, RNDINT, PLCI, NUMCI
      INTEGER    MAXCOL, COLEL, UNDF(3), EXPLCT, PERM(10), RTANS(10)
      INTEGER    NDISCT, NDLNST, PLCOL, RANLST

      REAL       THISLW, NOMLW,MINLW,MAXLW, YLOC,YINCR, XA(9),YA(9)
      PARAMETER (CSIZ = 6)
      REAL       BCKCOL(CSIZ), FORCOL(CSIZ), NEWCOL(CSIZ)
      REAL       TSTLW1,TSTLW2, MULT, RAD, CENTX,CENTY, PI

      CHARACTER  LDESCR(4)*14, PROMPT*200, LNTID*12
      CHARACTER  PFSW*1, MSG*2600, TSTYPE*1, SUFFIX*70

      INTEGER    IDUM1,IDUM2,IDUM3,IDUM4,IDUM5,IDUM6,IDUM7
      REAL       RDUM1,RDUM2,RDUM3

      LOGICAL    NOZERO

      DATA  LDESCR / 'solid', 'dashed', 'dotted', 'dashed-dotted' /

      IF (TSTYPE .EQ. 'i' .OR. TSTYPE .EQ. 'c') THEN
C        tstype valid
      ELSE
         CALL UNMSG ('PLATTR called with TSTYPE = ' // TSTYPE)
      ENDIF

C open PHIGS
      CALL XPOPPH (ERRFIL, MEMUN)
C set-up of workstation and dialogue area
      PICSTR = 101
      TXCI = 1
      CALL SETDLG (PICSTR, 801,TXCI)

      CALL POPST (PICSTR)
C  by convention, view #1 is for picture
      CALL PSVWI (1)
C  use individual attributes
      CALL SETASF (PINDIV)

C  use <inquire polyline facilities> to determine:
C    lavslt = list of available standard linetypes
C    lavrlt = list of available registered (non-mandatory) linetypes
C    lavilt = list of available implementor-defined linetypes
C             (may be derived directly from linetype value, if reported
C              number of linetypes < 0)
C    numlw  = number of available line-widths
C    nomlw  = nominal line-width (DC)
C    minlw,maxlw = minimum,maximum line-width (DC)
      CALL PQPLF (SPECWT,0, ERRIND, NUMLT,IDUM1,
     1            NUMLW,NOMLW,MINLW,MAXLW, IDUM2 )
      CALL CHKINQ ('pqplf', ERRIND)
      STDDX = 0
      REGDX = 0
      IMPDX = 0
      DO 100 LDX = 1,ABS(NUMLT)
         CALL PQPLF (SPECWT,LDX, ERRIND, IDUM1, THISLT,
     1               IDUM2,RDUM1,RDUM2,RDUM3, IDUM3)
         CALL CHKINQ ('pqplf', ERRIND)
         IF     (THISLT .GT. PLDASD) THEN
            REGDX = REGDX+1
            LAVRLT (REGDX) = THISLT
         ELSEIF (THISLT .GT. 0) THEN
            STDDX = STDDX+1
            LAVSLT (STDDX) = THISLT
         ELSE
            IF (IMPDX .LT. IMPMAX) THEN
               IF (NUMLT .LT. 0) THEN
C  build list from first LDX negative integers
                  IMPDX = IMPDX+1
                  LAVILT (IMPDX) = -LDX
               ELSE
                  IMPDX = IMPDX+1
                  LAVILT (IMPDX) = THISLT
               ENDIF
            ENDIF
         ENDIF
100   CONTINUE

      IF (TSTYPE .EQ. 'c') THEN
C  ndisct = number of distinct foreground colors in color table.
         CALL DISCOL (8, WKID, NDISCT)
      ENDIF

C  *** *** ***   1. linetype

C  mark start of linetypes
      CALL PLB (1)

C  Display in random order all linetypes in lavslt
      CALL RNPERM (STDDX, PERM)
      XA(1) = 0.3
      XA(2) = 0.9
      YINCR = 0.8/STDDX
      YLOC  = 0.9
      IF (TSTYPE .EQ. 'c') THEN
         CALL NDLW (YINCR/2, NOMLW,MINLW)
      ENDIF
      CALL NUMLAB (STDDX, 0.2, YLOC, YINCR)

      DO 200 LDX = 1,STDDX
         YA(1) = YLOC
         YA(2) = YLOC
         IF (TSTYPE .EQ. 'c') THEN
            CALL PSPLCI (RNDINT (1, NDISCT))
         ENDIF
         CALL PSLN (LAVSLT(PERM(LDX)))
         CALL PPL  (2, XA,YA)
         YLOC  = YLOC-YINCR
200   CONTINUE

C  mark end of linetypes
      CALL PLB (2)

      IF (TSTYPE .EQ. 'c') THEN
         SUFFIX = ', even when a non-default linewidth or color ' //
     1            'is used.'
      ELSE
         SUFFIX = '.'
      ENDIF

      CALL SETMSG ('3 4 5 6 8', 'The mandatory linetypes (1-4) ' //
     1             'should be recognizable from their standard ' //
     2             'description' // SUFFIX)

C  Operator must identify each linetype according to official
C  English description: solid, dotted, dashed, dotted-dashed.
      PROMPT = 'MANDATORY LINETYPES: List, in order, the numeric ' //
     1         'labels for linetypes: ' // LDESCR(LAVSLT(1))
      RTANS(PERM(1)) = 1
      DO 300 LDX = 2,STDDX
         PROMPT(ITRIM(PROMPT)+1:) = ', ' // LDESCR(LAVSLT(LDX))
         RTANS(PERM(LDX)) = LDX
300   CONTINUE
      PROMPT(ITRIM(PROMPT)+1:) = '.'

      CALL DLSTPF (PROMPT, STDDX, RTANS, 'L')

C  clear out last display from structure
      CALL PSEP (1)
      CALL PDELLB (1,2)

C look for unsupported negative linetype
      IF (NUMLT .GT. 0) THEN
C  find negative linetype *not* in LAVILT and whether zero
C  is supported.
         NOZERO = .TRUE.
         MINILT = 0
         DO 400 LDX = 1,IMPDX
            IF (LAVILT(LDX) .EQ. 0) THEN
               NOZERO = .FALSE.
            ENDIF
            IF (LAVILT(LDX) .LT. MINILT) THEN
               MINILT = LAVILT(LDX)
            ENDIF
400      CONTINUE
C  LUNSUP is list of unsupported linetypes, indexed by UNSUDX
         IF (NOZERO) THEN
            UNSUDX = 1
            LUNSUP (UNSUDX) = 0
         ELSE
            UNSUDX = 0
         ENDIF
C  since this value is less than minimum of LAVILT, it's unsupported
         UNSUDX = UNSUDX+1
         LUNSUP (UNSUDX) = MINILT - 1
      ELSE
C  all non-positive linetypes supported - derived directly from linetype value
         UNSUDX = 0
      ENDIF

C now look for unsupported positive linetype
      MAXILT = 5
      DO 500 LDX = 1,REGDX
         IF (LAVRLT(LDX) .GT. MAXILT) THEN
            MAXILT = LAVRLT(LDX)
         ENDIF
500   CONTINUE
      UNSUDX = UNSUDX+1
      LUNSUP (UNSUDX) = MAXILT + 1

C  now add to end of LUNSUP entry #1 and a non-solid linetype
      UNSUDX = UNSUDX+1
      LUNSUP (UNSUDX) = 1

C  other = another available non-solid linetype: 2,3, or 4
      DO 600 LDX = 1,STDDX
         IF (LAVSLT(LDX) .GT. 1) THEN
C  this is the only one that should appear as non-solid
            UNSUDX = UNSUDX+1
            LUNSUP (UNSUDX) = LAVSLT(LDX)
            GOTO 610
         ENDIF
600   CONTINUE
C  no other available  - goto end_solid
      GOTO 690

610   CONTINUE

C  Display in random order linetypes in LUNSUP
      CALL RNPERM (UNSUDX, PERM)
      XA(1) = 0.3
      XA(2) = 0.9
      YINCR = 0.8/UNSUDX
      YLOC  = 0.9
      IF (TSTYPE .EQ. 'c') THEN
         CALL NDLW (YINCR/2, NOMLW,MINLW)
      ENDIF
      CALL NUMLAB (UNSUDX, 0.2, YLOC,YINCR)

      DO 700 LDX = 1,UNSUDX
         YA(1) = YLOC
         YA(2) = YLOC
         IF (TSTYPE .EQ. 'c') THEN
            CALL PSPLCI (RNDINT (1, NDISCT))
         ENDIF
         CALL PSLN (LUNSUP(PERM(LDX)))
         CALL PPL  (2, XA,YA)
         YLOC  = YLOC-YINCR
700   CONTINUE

      CALL SETMSG ('3 4 5 6 7 8', 'Unavailable linetypes should be ' //
     1             'displayed as linetype number 1' // SUFFIX)
      CALL DCHPF ('UNSUPPORTED LINETYPES: Which line is NOT solid?',
     1             UNSUDX, UNSUDX, PERM)

C  clear out last display from structure
      CALL PSEP (1)
      CALL PDELLB (1,2)

C  end_solid:
690   CONTINUE

C  REGDX = size of lavrlt
      IF (REGDX .LT. 1) THEN
         CALL INMSG ('No non-mandatory registered linetypes are ' //
     1               'supported.')
         GOTO 890
      ENDIF
C display and label at most 8 linetypes in lavrlt
      NUMDIS = MIN(8, REGDX)
      CALL RNSET (NUMDIS, REGDX, PERM)

C  build list of linetypes to be displayed
      DO 770 LDX = 1,NUMDIS
         LINDIS(LDX) = LAVRLT(PERM(LDX))
770   CONTINUE
C  and sort it
      CALL SRTIAR (NUMDIS, LINDIS)

C  Display in order (some) linetypes in LAVRLT
      XA(1) = 0.5
      XA(2) = 0.9
      YINCR = 0.8/NUMDIS
      YLOC  = 0.9
      IF (TSTYPE .EQ. 'c') THEN
         CALL NDLW (YINCR/2, NOMLW,MINLW)
      ENDIF

      DO 800 LDX = 1,NUMDIS
         YA(1) = YLOC
         YA(2) = YLOC
         THISLT = LINDIS(LDX)
         WRITE (LNTID, '(A,I2,A)') 'Type#', THISLT, ':'
         CALL VISLAB (LNTID(1:8), 'R', 0.05,0.45,
     1                YLOC-0.3*YINCR, YLOC+0.3*YINCR)
         IF (TSTYPE .EQ. 'c') THEN
            CALL PSPLCI (RNDINT (1, NDISCT))
         ENDIF
         CALL PSLN (THISLT)
         CALL PPL  (2, XA,YA)
         YLOC  = YLOC-YINCR
800   CONTINUE
      CALL SETMSG ('3 4 5 6 9', 'The registered non-mandatory '       //
     1             'linetypes (> 4) should agree with their '         //
     2             'registered description' // SUFFIX)
      CALL DYNPF ('NON-MANDATORY REGISTERED LINETYPES: Is each '  //
     1            'linetype depicted according to its numeric '   //
     2            'identifier''s specification in the ISO register?',
     3            'Y')

C  clear out last display from structure
      CALL PSEP (1)
      CALL PDELLB (1,2)

C  neg_type:
890   CONTINUE

C  IMPDX = size of lavilt
      IF (IMPDX .LT. 1) THEN
         CALL INMSG ('No implementor-defined linetypes are supported.')
         GOTO 990
      ENDIF
C display and label at most 8 linetypes in lavilt
      NUMDIS = MIN(8, IMPDX)
      CALL RNSET (NUMDIS, IMPDX, PERM)

C  build list of linetypes to be displayed
      DO 870 LDX = 1,NUMDIS
         LINDIS(LDX) = LAVILT(PERM(LDX))
870   CONTINUE
C  and sort it
      CALL SRTIAR (NUMDIS, LINDIS)

C  Display (some) linetypes in LAVILT
      XA(1) = 0.5
      XA(2) = 0.9
      YINCR = 0.8/NUMDIS
      YLOC  = 0.9
      IF (TSTYPE .EQ. 'c') THEN
         CALL NDLW (YINCR/2, NOMLW,MINLW)
      ENDIF

      DO 900 LDX = 1,NUMDIS
         YA(1) = YLOC
         YA(2) = YLOC
         THISLT = LINDIS(LDX)
         WRITE (LNTID, '(A,I5,A)') 'Type#', THISLT, ':'
         CALL VISLAB (LNTID(1:11), 'R', 0.05,0.45,
     1                YLOC-0.3*YINCR, YLOC+0.3*YINCR)
         IF (TSTYPE .EQ. 'c') THEN
            CALL PSPLCI (RNDINT (1, NDISCT))
         ENDIF
         CALL PSLN (THISLT)
         CALL PPL  (2, XA,YA)
         YLOC  = YLOC-YINCR
900   CONTINUE
      CALL SETMSG ('3 4 5 6 10', 'Implementor-defined linetypes (< ' //
     1             '1) should agree with the descriptions supplied ' //
     2             'by the implementor' // SUFFIX)
      CALL DYNPF ('IMPLEMENTOR DEFINED LINETYPES: Is each linetype ' //
     1            'depicted according to its numeric identifier''s ' //
     2            'specification in the implementor documentation?',
     3            'Y')

C  clear out last display from structure
      CALL PSEP (1)
      CALL PDELLB (1,2)

C  Check if all implementor linetypes already shown or
C  linetype appearance derived directly from value - if so, skip test
      IF (IMPDX .LE. 8 .OR. NUMLT .LT. 0) GOTO 990
      IF (TSTYPE .EQ. 'c') GOTO 990

      CALL SETMSG ('10', 'All reported implementor-defined '      //
     1             'linetypes available for a given workstation ' //
     2             'type should be documented by the implementor.')

      CALL SRTIAR (IMPDX, LAVILT)
      WRITE (MSG, '(500I5)') (LAVILT(LDX), LDX=1,IMPDX)
      CALL DYNPF ('COMPLETE DOCUMENTATION FOR IMPLEMENTOR LINETYPES: '//
     1            'Are all these reportedly available linetypes '     //
     2            'documented by the implementor?' // MSG, 'Y')

C  end_linetype:
990   CONTINUE

C  *** *** ***   2. line-width

      IF (TSTYPE .EQ. 'c') THEN
         SUFFIX = ', even when a non-default linetype ' //
     3            'or color is used.'
      ELSE
         SUFFIX = '.'
         NDLNST = 1
         PLCOL  = 1
      ENDIF

C  Are there at least two visually distinguishable line-widths?
      IF (NUMLW .EQ. 1          .OR.
     1    MAXLW .LE. 1.02*MINLW .OR.
     2    MAXLW-MINLW .LT. QVIS)    THEN
         TSTLW1 = MAXLW
         TSTLW2 = MAXLW
         MULT   = 2.0
      ELSE
         TSTLW1 = MINLW
         TSTLW2 = MAXLW
         IF (NUMLW .EQ. 0) THEN
C  continuous range of linewidths available - take 4 geometric steps
            MULT = 0.9999 * ((TSTLW2/TSTLW1) ** 0.25)
         ELSE
C  discrete set of linewidths available - test min/max only
            MULT = 0.9999 * (TSTLW2/TSTLW1)
         ENDIF
      ENDIF

C test scale sizes
C set up counter for number of sizes passed
      NUMPAS = 0
      THISLW = TSTLW1
C  loop thru various linewidths/next_lw:
1000  CONTINUE

      IF (TSTYPE .EQ. 'c') THEN
C  random linetype and color:
         NDLNST = RANLST (STDDX,LAVSLT, REGDX,LAVRLT, IMPDX,LAVILT)
         PLCOL  = RNDINT (1, NDISCT)
      ENDIF
C  invoke SHOWLW subroutine to test requested line width = thislw
C                                    expected line width = thislw
      CALL SHOWLW ('VARIOUS LINEWIDTH SCALE FACTORS', NDLNST, PLCOL,
     1              THISLW, THISLW, NOMLW, MINLW, PFSW)

      IF (PFSW .EQ. 'A') THEN
         IF (NUMPAS .LT. 2) THEN
            CALL INMSG ('Skipping test case for various scale '     //
     1                  'factors because largest valid line width ' //
     2                  'is greater than screen size.')
            GOTO 1110
         ELSE
            GOTO 1100
         ENDIF
      ELSEIF (PFSW .NE. 'P') THEN
         NUMPAS = 0
         GOTO 1100
      ENDIF
C  if PFSW = 'P', set up for next linewidth to be tested:
      NUMPAS = NUMPAS+1
      THISLW = THISLW * MULT
      IF (THISLW .LE. TSTLW2) GOTO 1000

1100  CONTINUE
      CALL SETMSG ('3 13 14 15', 'Available linewidth scale '       //
     1             'factors should control the realized thickness ' //
     2             'of a polyline' // SUFFIX)
      CALL IFPF (NUMPAS.GT.0)

C  min_max_coerce:
1110  CONTINUE
      CALL SETMSG ('3 13 14 16', 'A requested positive linewidth '    //
     1             'scale factor below the minimum available should ' //
     2             'be realized as the minimum' // SUFFIX)

      IF (TSTYPE .EQ. 'c') THEN
C  random linetype and color:
         NDLNST = RANLST (STDDX,LAVSLT, REGDX,LAVRLT, IMPDX,LAVILT)
         PLCOL  = RNDINT (1, NDISCT)
      ENDIF
      CALL SHOWLW ('POSITIVE LINEWIDTH BELOW MINIMUM', NDLNST, PLCOL,
     1              MINLW/2, MINLW, NOMLW, MINLW, PFSW)
      CALL IFPF (PFSW .EQ. 'P')

      CALL SETMSG ('3 13 14 16', 'A requested negative linewidth '  //
     1             'scale factor should be realized as the minimum' //
     2             SUFFIX)
      IF (TSTYPE .EQ. 'c') THEN
C  random linetype and color:
         NDLNST = RANLST (STDDX,LAVSLT, REGDX,LAVRLT, IMPDX,LAVILT)
         PLCOL  = RNDINT (1, NDISCT)
      ENDIF
      CALL SHOWLW ('NEGATIVE LINEWIDTH', NDLNST, PLCOL,
     1             -MAXLW-100, MINLW, NOMLW, MINLW, PFSW)
      CALL IFPF (PFSW .EQ. 'P')

C test line width scale above maximum = maximum
      IF (TSTYPE .EQ. 'c') THEN
C  random linetype and color:
         NDLNST = RANLST (STDDX,LAVSLT, REGDX,LAVRLT, IMPDX,LAVILT)
         PLCOL  = RNDINT (1, NDISCT)
      ENDIF
      CALL SHOWLW ('POSITIVE LINEWIDTH ABOVE MAXIMUM', NDLNST, PLCOL,
     1              MAXLW*2, MAXLW, NOMLW, MINLW, PFSW)
      IF (PFSW .EQ. 'A') THEN
         CALL INMSG ('Skipping test case for requested scale factor ' //
     1               'greater than maximum because largest valid '    //
     2               'line width is greater than screen size.')
      ELSE
         CALL SETMSG ('3 13 14 16', 'A requested linewidth scale '    //
     1                'factor above the maximum available should be ' //
     2                'realized as the maximum' // SUFFIX)
         CALL IFPF (PFSW .EQ. 'P')
      ENDIF

C  *** *** ***   3. polyline color index

      IF (TSTYPE .EQ. 'c') THEN
         SUFFIX = ', even when a non-default linetype ' //
     3            'or linewidth is used.'
      ELSE
         SUFFIX = '.'
      ENDIF

      CALL SETMSG ('3 19 20 21', 'A defined polyline color index '  //
     1             'should cause the addressed entry in the color ' //
     2             'table to be used when rendering a polyline'     //
     3             SUFFIX)

      CALL PSCMD (WKID, PRGB)
C  bckcol = background color = realized color spec for entry #0
C  forcol = foreground color = realized color spec for entry #1
      CALL PQCR (WKID, 0, CSIZ, PREALI, ERRIND, CDIM, BCKCOL)
      CALL CHKINQ ('pqcr', ERRIND)
      CALL PQCR (WKID, 1, CSIZ, PREALI, ERRIND, CDIM, FORCOL)
      CALL CHKINQ ('pqcr', ERRIND)

C  szcolt = maximum size of color table (including entry #0)
      CALL PQWKSL (SPECWT, ERRIND, IDUM1,IDUM2,IDUM3,IDUM4,IDUM5,IDUM6,
     1             SZCOLT, IDUM7)
      CALL CHKINQ ('pqwksl', ERRIND)
C  how many lines to draw
      NUMLIN = MIN(8, SZCOLT)
C  from values 2 to SZCOLT-1, pick NUMLIN-2 entries
      CALL RNBSET (NUMLIN-2, 2, SZCOLT-1, LNCOL)
C  add in 0 and 1
      LNCOL(NUMLIN-1) = 1
      LNCOL(NUMLIN)   = 0
      CALL RNPERM (NUMLIN, PERM)

C  from lncol, randomly select entry #visdx, but not the one that
C    contains 0, since this may not be re-settable.
      VISDX  = RNDINT(1,NUMLIN-1)
      VISCOL = LNCOL(VISDX)

C  ensure dialogue visible
      DTXCI  = VISCOL
      CALL PSTXCI (VISCOL)
      YLOC  = 0.9
      YINCR = 0.8 / (NUMLIN-1)
      XA(1) = .2
      XA(2) = .9

      IF (TSTYPE .EQ. 'c') THEN
C  non-default linewidth
         CALL NDLW (YINCR/2, NOMLW,MINLW)
      ENDIF
      CALL NUMLAB (NUMLIN, .15, YLOC, YINCR)

C  for each color entry - set to background color and draw line
      DO 2100 IX = 1,NUMLIN
         PLCI = LNCOL(PERM(IX))
         CALL PSCR (WKID, PLCI, CDIM, BCKCOL)
         CALL PSPLCI (PLCI)
         IF (TSTYPE .EQ. 'c') THEN
C  random linetype
            CALL PSLN (RANLST(STDDX,LAVSLT, REGDX,LAVRLT, IMPDX,LAVILT))
         ENDIF
         YA(1) = YLOC
         YA(2) = YLOC
         CALL PPL (2, XA,YA)
         YLOC  = YLOC-YINCR
2100  CONTINUE

C  set entry lncol(visdx) to forcol
      CALL PSCR (WKID, VISCOL, CDIM, FORCOL)

      CALL DCHPF ('DEFINED POLYLINE COLOR INDICES: Which line ' //
     1            'is visible?', NUMLIN, VISDX, PERM)

C  clear out last display from structure
      CALL PSEP (1)
      CALL PDELLB (1,2)

      CALL SETMSG ('3 19 20 22 ', 'An undefined polyline color '     //
     1             'index should cause entry number 1 in the color ' //
     2             'table to be used when rendering a polyline'      //
     3             SUFFIX)

C  Try to set entry #1 opposite from BCKCOL
      DO 2150 IX = 1,3
         IF (BCKCOL(IX) .GT. 0.5) THEN
            NEWCOL(IX) = 0.0
         ELSE
            NEWCOL(IX) = 1.0
         ENDIF
2150  CONTINUE

C  set entry #1 different from FORCOL - make sure other default
C  is to current color-rep of #1, not just a predefined color.
      IF (FORCOL(3) .GT. 0.5) THEN
         NEWCOL(3) = 0.0
      ELSE
         NEWCOL(3) = 1.0
      ENDIF

      CALL PSCR  (WKID, 1, CDIM, NEWCOL)

      CALL PQECI (WKID, 0, ERRIND, NUMCI, IDUM1)
      CALL CHKINQ ('pqeci', ERRIND)
      MAXCOL = 0
      DO 2200 IX = 1,NUMCI
         CALL PQECI (WKID, IX, ERRIND, IDUM1, COLEL)
         CALL CHKINQ ('pqeci', ERRIND)
         MAXCOL = MAX(MAXCOL, COLEL)
2200  CONTINUE
C  undf1,undf2,undf3 = 3 undefined, positive color indices
      UNDF(1) = MAXCOL+1
      UNDF(2) = MAXCOL+10
      UNDF(3) = MAXCOL+90

C  explct = number of explicit lines of color #1 = random integer
C  from 0 to 4
      EXPLCT = RNDINT(0,4)
      NUMLIN = EXPLCT+3
      CALL RNPERM (NUMLIN, PERM)

C  draw star with color #1
      RAD = .15
      CENTX = .5
      CENTY = .75
      PI = 3.14159265
      DO 2400 IX = 1,6
         YA(IX) = CENTY + RAD*COS((4*PI*IX)/5)
         XA(IX) = CENTX + RAD*SIN((4*PI*IX)/5)
2400  CONTINUE
      CALL PSPLCI (1)
      CALL PPL (6,XA,YA)

C  display interleaved:
C    three lines of color u1,u2,u3,
C    explct lines of color #1
      YLOC  = 0.5
      YINCR = 0.5 / 8
      XA(1) = .1
      XA(2) = .9

      IF (TSTYPE .EQ. 'c') THEN
C  non-default linewidth
         CALL NDLW (YINCR/2, NOMLW,MINLW)
      ENDIF

C  for each color entry - set to background color and draw line
      DO 2300 IX = 1,NUMLIN
         COLEL = PERM(IX)
         IF (COLEL .LE. 3) THEN
            CALL PSPLCI (UNDF(COLEL))
         ELSE
            CALL PSPLCI (1)
         ENDIF
         IF (TSTYPE .EQ. 'c') THEN
C  random linetype
            CALL PSLN (RANLST(STDDX,LAVSLT, REGDX,LAVRLT, IMPDX,LAVILT))
         ENDIF
         YA(1) = YLOC
         YA(2) = YLOC
         CALL PPL (2, XA,YA)
         YLOC  = YLOC-YINCR
2300  CONTINUE

      CALL DCHPFV ('UNDEFINED POLYLINE COLOR INDICES: How many '    //
     1             'of the horizontal lines are the same color as ' //
     2             'the star?', 12, NUMLIN)

666   CONTINUE
C  wrap it up.
      END


text:
      SUBROUTINE PZLACPY( UPLO, M, N, A, IA, JA, DESCA, B, IB, JB,
     $                    DESCB )
*
*  -- ScaLAPACK auxiliary routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     May 1, 1997
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            IA, IB, JA, JB, M, N
*     ..
*     .. Array Arguments ..
      INTEGER            DESCA( * ), DESCB( * )
      COMPLEX*16         A( * ), B( * )
*     ..
*
*  Purpose
*  =======
*
*  PZLACPY copies all or part of a distributed matrix A to another
*  distributed matrix B.  No communication is performed, PZLACPY
*  performs a local copy sub( A ) := sub( B ), where sub( A ) denotes
*  A(IA:IA+M-1,JA:JA+N-1) and sub( B ) denotes B(IB:IB+M-1,JB:JB+N-1).
*
*  Notes
*  =====
*
*  Each global data object is described by an associated description
*  vector.  This vector stores the information required to establish
*  the mapping between an object element and its corresponding process
*  and memory location.
*
*  Let A be a generic term for any 2D block cyclicly distributed array.
*  Such a global array has an associated description vector DESCA.
*  In the following comments, the character _ should be read as
*  "of the global array".
*
*  NOTATION        STORED IN      EXPLANATION
*  --------------- -------------- --------------------------------------
*  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
*                                 DTYPE_A = 1.
*  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
*                                 the BLACS process grid A is distribu-
*                                 ted over. The context itself is glo-
*                                 bal, but the handle (the integer
*                                 value) may vary.
*  M_A    (global) DESCA( M_ )    The number of rows in the global
*                                 array A.
*  N_A    (global) DESCA( N_ )    The number of columns in the global
*                                 array A.
*  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
*                                 the rows of the array.
*  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
*                                 the columns of the array.
*  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
*                                 row of the array A is distributed.
*  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
*                                 first column of the array A is
*                                 distributed.
*  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
*                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
*
*  Let K be the number of rows or columns of a distributed matrix,
*  and assume that its process grid has dimension p x q.
*  LOCr( K ) denotes the number of elements of K that a process
*  would receive if K were distributed over the p processes of its
*  process column.
*  Similarly, LOCc( K ) denotes the number of elements of K that a
*  process would receive if K were distributed over the q processes of
*  its process row.
*  The values of LOCr() and LOCc() may be determined via a call to the
*  ScaLAPACK tool function, NUMROC:
*          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
*          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
*  An upper bound for these quantities may be computed by:
*          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
*          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
*
*  Arguments
*  =========
*
*  UPLO    (global input) CHARACTER
*          Specifies the part of the distributed matrix sub( A ) to be
*          copied:
*          = 'U':   Upper triangular part is copied; the strictly
*                   lower triangular part of sub( A ) is not referenced;
*          = 'L':   Lower triangular part is copied; the strictly
*                   upper triangular part of sub( A ) is not referenced;
*          Otherwise:  All of the matrix sub( A ) is copied.
*
*  M       (global input) INTEGER
*          The number of rows to be operated on i.e the number of rows
*          of the distributed submatrix sub( A ). M >= 0.
*
*  N       (global input) INTEGER
*          The number of columns to be operated on i.e the number of
*          columns of the distributed submatrix sub( A ). N >= 0.
*
*  A       (local input) COMPLEX*16 pointer into the local memory
*          to an array of dimension (LLD_A, LOCc(JA+N-1) ). This array
*          contains the local pieces of the distributed matrix sub( A )
*          to be copied from.
*
*  IA      (global input) INTEGER
*          The row index in the global array A indicating the first
*          row of sub( A ).
*
*  JA      (global input) INTEGER
*          The column index in the global array A indicating the
*          first column of sub( A ).
*
*  DESCA   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix A.
*
*  B       (local output) COMPLEX*16 pointer into the local memory
*          to an array of dimension (LLD_B, LOCc(JB+N-1) ). This array
*          contains on exit the local pieces of the distributed matrix
*          sub( B ) set as follows:
*
*          if UPLO = 'U', B(IB+i-1,JB+j-1) = A(IA+i-1,JA+j-1),
*                         1<=i<=j, 1<=j<=N;
*          if UPLO = 'L', B(IB+i-1,JB+j-1) = A(IA+i-1,JA+j-1),
*                         j<=i<=M, 1<=j<=N;
*          otherwise,     B(IB+i-1,JB+j-1) = A(IA+i-1,JA+j-1),
*                         1<=i<=M, 1<=j<=N.
*
*  IB      (global input) INTEGER
*          The row index in the global array B indicating the first
*          row of sub( B ).
*
*  JB      (global input) INTEGER
*          The column index in the global array B indicating the
*          first column of sub( B ).
*
*  DESCB   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix B.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
     $                   LLD_, MB_, M_, NB_, N_, RSRC_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
     $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
     $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, IAA, IBB, IBLK, IN, ITMP, J, JAA, JBB,
     $                   JBLK, JN, JTMP
*     ..
*     .. External Subroutines ..
      EXTERNAL           PZLACP2
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ICEIL
      EXTERNAL           ICEIL, LSAME
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MIN, MOD
*     ..
*     .. Executable Statements ..
*
      IF( M.EQ.0 .OR. N.EQ.0 )
     $   RETURN
*
      IN = MIN( ICEIL( IA, DESCA( MB_ ) ) * DESCA( MB_ ), IA+M-1 )
      JN = MIN( ICEIL( JA, DESCA( NB_ ) ) * DESCA( NB_ ), JA+N-1 )
*
      IF( M.LE.( DESCA( MB_ ) - MOD( IA-1, DESCA( MB_ ) ) ) .OR.
     $    N.LE.( DESCA( NB_ ) - MOD( JA-1, DESCA( NB_ ) ) ) ) THEN
         CALL PZLACP2( UPLO, M, N, A, IA, JA, DESCA,
     $                 B, IB, JB, DESCB )
      ELSE
*
         IF( LSAME( UPLO, 'U' ) ) THEN
            CALL PZLACP2( UPLO, IN-IA+1, N, A, IA, JA, DESCA,
     $                    B, IB, JB, DESCB )
            DO 10 I = IN+1, IA+M-1, DESCA( MB_ )
               ITMP = I-IA
               IBLK = MIN( DESCA( MB_ ), M-ITMP )
               IBB = IB + ITMP
               JBB = JB + ITMP
               JAA = JA + ITMP
               CALL PZLACP2( UPLO, IBLK, N-ITMP, A, I, JAA, DESCA,
     $                       B, IBB, JBB, DESCB )
   10       CONTINUE
         ELSE IF( LSAME( UPLO, 'L' ) ) THEN
            CALL PZLACP2( UPLO, M, JN-JA+1, A, IA, JA, DESCA,
     $                    B, IB, JB, DESCB )
            DO 20 J = JN+1, JA+N-1, DESCA( NB_ )
               JTMP = J-JA
               JBLK = MIN( DESCA( NB_ ), N-JTMP )
               IBB = IB + JTMP
               JBB = JB + JTMP
               IAA = IA + JTMP
               CALL PZLACP2( UPLO, M-JTMP, JBLK, A, IAA, J, DESCA,
     $                       B, IBB, JBB, DESCB )
   20       CONTINUE
         ELSE
            IF( M.LE.N ) THEN
               CALL PZLACP2( UPLO, IN-IA+1, N, A, IA, JA, DESCA,
     $                       B, IB, JB, DESCB )
               DO 30 I = IN+1, IA+M-1, DESCA( MB_ )
                  ITMP = I-IA
                  IBLK = MIN( DESCA( MB_ ), M-ITMP )
                  IBB = IB+ITMP
                  CALL PZLACP2( UPLO, IBLK, N, A, I, JA, DESCA,
     $                          B, IBB, JB, DESCB )
   30          CONTINUE
            ELSE
               CALL PZLACP2( UPLO, M, JN-JA+1, A, IA, JA, DESCA,
     $                       B, IB, JB, DESCB )
               DO 40 J = JN+1, JA+N-1, DESCA( NB_ )
                  JTMP = J-JA
                  JBLK = MIN( DESCA( NB_ ), N-JTMP )
                  JBB = JB+JTMP
                  CALL PZLACP2( UPLO, M, JBLK, A, IA, J, DESCA,
     $                          B, IB, JBB, DESCB )
   40          CONTINUE
            END IF
         END IF
*
      END IF
*
      RETURN
*
*     End of PZLACPY
*
      END

ext: f

max_stars_repo_name: jschueller/scalapack

max_issues_repo_name: jschueller/scalapack

max_forks_repo_name: jschueller/scalapack


text:
!****************************************************************************
!
!  PROGRAM: web
!
!  PURPOSE:  Main shell program.
!
!****************************************************************************

      program versatile_predict
      use file_io
      use constituent
	  implicit none
	  integer , parameter :: KIN = 4
      character(LEN=80) :: main_file
      integer, parameter :: max_main_cnstnt = 80
      integer, parameter :: max_infer = 80
      integer :: ndef
      integer :: nloc
      integer :: itrend
      integer :: id1,im1,iy1,id2,im2,iy2,ic1,ic2

      integer, allocatable :: ipredtime(:,:)

!c    For prediction
      character*80 :: file_predict
      integer :: station_no
      integer kd1,kd2
      character*20 :: station_name_pred
      integer :: latd,latm
      integer :: lond,lonm
      integer :: nconstituent
      integer :: ntvl_min
      character*5, dimension(max_main_cnstnt) :: all_names !todo: should be all constituents
      real*8, allocatable,dimension(:) :: all_freq
      real*8, allocatable,dimension(:,:,:) :: all_amp
      real*8, allocatable,dimension(:,:,:) :: all_phase
      real*8,allocatable :: values(:,:,:)
      integer i
      integer :: idd,imm,iyy,icc,ikd
      integer nperday,istep,imin,ihh,isec,ntime
      character :: sep1,sep2
      integer ::  istart_ha(7)  ! timestamp of original analysis start (for trend)
      integer ::  iend_ha(7)    ! timestamp of original analysis end

!c	FILE I/O
!c	KIN is the master input file.
!c	fnam6 is the file to which the output is sent. It is assigned the number
!c		lp.
!c	fnam7 is file containing the constituents to be included in the analysis,
!c		the analysis period, inference parameters, the flag controlling
!c		height or current analyses, and site information. It is assigned the
!c		number kr1.
!c	fnam8 is the file containing all the astronomical argument information
!c		(it should not have to be changed)
!c	fnam9 is the file containing the observations and their times. It is
!c		assigned the number kr2.
!c	fnam11 is a file to which information on the SVD matrix fit is output when
!c		ibin > 0. This info includes matrix covariances which are useful
!c		in determining (in)dependence of the chosen constituents.
!c	Original, fitted, and residual time series are output to file 25 while
!c	the same are also output to file 26 in a format that could be input to
!c	Excel for plotting.
!c
      call parse_command_line(main_file)

      open(UNIT=KIN,FILE=main_file,STATUS='OLD',ACTION='READ')
      read(KIN,*) file_constituent_out
      read(KIN,*) file_const_db
      read(KIN,*) file_predict
      read(KIN,'(I2,I2,1a,I2,1a,I2)') ic1,iy1,sep1,im1,sep2,id1
      read(KIN,'(I2,I2,1a,I2,1a,I2)') ic2,iy2,sep1,im2,sep2,id2
      read(KIN,*) ntvl_min
      close(KIN)

      call GDAY(ID1,IM1,IY1,IC1,kd1)
      call GDAY(ID2,IM2,IY2,IC2,kd2)
      nperday=(24*60)/ntvl_min          ! todo: mention noneven intervals
      ntime = ((kd2-kd1) + 1) * nperday ! todo: is this calculation correct?
      if (.not. allocated(ipredtime)) allocate(ipredtime(7,ntime))

      i=0
      do ikd=kd1,kd2
          call dmy(IDD,IMM,IYY,ICC,ikd)
          do istep=1,nperday
              i=i+1
              imin=(istep-1)*ntvl_min
              ihh=imin/60
              imin=mod(imin,60)
              isec = 0
              ipredtime(1,i)=iyy
              ipredtime(2,i)=imm
              ipredtime(3,i)=idd
              ipredtime(4,i)=ihh
              ipredtime(5,i)=imin
              ipredtime(6,i)=isec
	        ipredtime(7,i)=icc
          end do
      end do
      ntime=i

      call vt_read_constituent_db(file_const_db)

      ! At the moment we are really only calling this to get ndef and nloc for allocation
      ! todo: could make the number of constituents dynamic as well.
      call vt_read_constituent_file_header(file_constituent_out,
     1                         station_no,
     1                         station_name_pred,
     1                         latd, latm,
     1                         lond,lonm,
     1                         istart_ha,iend_ha,
     1                         nconstituent,ndef,nloc,itrend)

      if (.not. allocated(all_freq)) allocate(all_freq(max_main_cnstnt))
      if (.not. allocated(all_amp)) allocate(all_amp(max_main_cnstnt,ndef,nloc))
      if (.not. allocated(all_phase)) allocate(all_phase(max_main_cnstnt,ndef,nloc))

      call vt_read_constituent_file(file_constituent_out,
     1                         station_no,
     1                         station_name_pred,
     1                         latd, latm,
     1                         lond,lonm,
     1                         istart_ha,
     1                         iend_ha,
     1                         nconstituent,ndef,nloc,itrend,
     1                         all_names,
     1                         all_freq,
     1                         all_amp,
     1                         all_phase)

      if (.not. allocated(values)) allocate(values(ntime,ndef,nloc))

      call vt_predict(ipredtime,
     1                ntime,
     1                values,
     1                ndef,
     1                nloc,
     1                latd, latm,
     1                istart_ha,
     1                iend_ha,
     1                nconstituent,
     1                itrend,
     1                all_names,
     1                all_freq,
     1                all_amp,
     1                all_phase)

      call write_fit_time_series(file_predict,ipredtime,values,ntime,ndef,nloc)

      deallocate(all_freq)
      deallocate(all_phase)
      deallocate(all_amp)
 !c     deallocate(obsdata)
	end program versatile_predict

      subroutine parse_command_line(buffer)
      character(LEN=*) buffer
      integer status
      call get_command_argument(1,buffer,status)
      buffer = trim(buffer)
      return
      end subroutine

ext: for

max_stars_repo_name: kjnam/vtide

max_issues_repo_name: esatel/vtide

max_forks_repo_name: esatel/vtide


text:
      SUBROUTINE PHICOR(I,DPHI1,DPHI2)
*
*
*       Differential potential corrections.
*       -----------------------------------
*
      INCLUDE 'common6.h'
*
*
*       Distinguish between single particle and (un)perturbed binary c.m.
      DPHI1 = 0.0
      DPHI2 = 0.0
      IF (I.LE.N) THEN
          NNB = LIST(1,I)
*       Search neighbour list backwards until the last c.m. particle.
          DO 10 L = NNB+1,2,-1
              J = LIST(L,I)
              IF (J.LE.N) GO TO 40
              RIJ2 = 0.0
              DO 5 K = 1,3
                  RIJ2 = RIJ2 + (X(K,I) - X(K,J))**2
    5         CONTINUE
*       Perform differential correction inside standard c.m. separation.
              IF (RIJ2.LT.CMSEP2*RMIN2) THEN
                  RJ1 = 0.0
                  RJ2 = 0.0
                  J1 = 2*(J - N) - 1
                  J2 = J1 + 1
                  DO 8 K = 1,3
                      RJ1 = RJ1 + (X(K,I) - X(K,J1))**2
                      RJ2 = RJ2 + (X(K,I) - X(K,J2))**2
    8             CONTINUE
                  DPHI1 = DPHI1 - BODY(J)/SQRT(RIJ2) +
     &                           BODY(J1)/SQRT(RJ1) + BODY(J2)/SQRT(RJ2)
              END IF
   10     CONTINUE
      ELSE
          I1 = 2*(I - N) - 1
          I2 = I1 + 1
          NP = LIST(1,I1)
*       Loop over all perturbers (only c.m. approximation for singles).
          DO 30 L = 2,NP+1
              J = LIST(L,I1)
              RIJ2 = 0.0
              DO 12 K = 1,3
                  RIJ2 = RIJ2 + (X(K,I) - X(K,J))**2
   12         CONTINUE
*       Treat case of close perturbing binary more carefully below.
              IF (J.LE.N.OR.RIJ2.GT.CMSEP2*RMIN2) THEN
                  RI1 = 0.0
                  RI2 = 0.0
                  DO 15 K = 1,3
                      RI1 = RI1 + (X(K,I1) - X(K,J))**2
                      RI2 = RI2 + (X(K,I2) - X(K,J))**2
   15             CONTINUE
*       Subtract c.m. interaction for each component (copied from GPU).
                  DPHI1 = DPHI1 - BODY(J)/SQRT(RIJ2) + BODY(J)/SQRT(RI1)
                  DPHI2 = DPHI2 - BODY(J)/SQRT(RIJ2) + BODY(J)/SQRT(RI2)
              ELSE
                  DPHI1 = DPHI1 - BODY(J)/SQRT(RIJ2)
                  DPHI2 = DPHI2 - BODY(J)/SQRT(RIJ2)
*       Include pair-wise interactions after subtracting the c.m.-c.m. terms.
                  J1 = 2*(J - N) - 1
                  J2 = J1 + 1
   25             RJ1 = 0.0
                  RJ2 = 0.0
                  DO 28 K = 1,3
                      RJ1 = RJ1 + (X(K,I1) - X(K,J1))**2
                      RJ2 = RJ2 + (X(K,I2) - X(K,J1))**2
   28             CONTINUE
*       Add each interaction separately (J1 followed by J2).
                  DPHI1 = DPHI1 + BODY(J1)/SQRT(RJ1)
                  DPHI2 = DPHI2 + BODY(J1)/SQRT(RJ2)
                  J1 = J1 + 1
                  IF (J1.EQ.J2) GO TO 25
              END IF
   30     CONTINUE
      END IF
*
   40 RETURN
*
      END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
Curves is a fitness circuit training program designed for busy women teenaged and older.  The local Curves is part of a large gyms gym chain, Curves has over one franchise location per every two McDonalds in the United States. There are 10,000 locations in 67 countries to work out in.  A complete workout takes just 30 minutes because of the efficiency of the program. They base their membership in part on gender.  Their intent is to create a welcoming environment for women.  They are also famed for their no mirrors, no makeup policy.

Free Weight Management and Nutrition Classes are also available.  Curves Circuit with Zumba has been available since January, 2011.  It is a perfect combination of strength training for all muscles and bones on the Curves hydraulic resistance equipment with the exciting fun cardio of Zumba.

The local Curves fitness center is a franchise. This local fitness center is owned by a local owner who pays the Curves company to use the Curves name. The Curves company gets a percentage of the profits.

 Community Service

 Curves donated hundreds of dollars of personal items for Mothers Day Baskets of Love project in May. This is the 4th year of supporting and encouraging the women who are in the homeless and domestic violence shelters of Yolo County.  32 baskets were filled with chocolate, tea & coffee, toiletries, jewelry, journals, books, mugs, makeup, and gifts to brighten their lives on this special day.

 Curves donated over 3,000 pounds of food and dollars to STEAC, our local food bank, in the 14th annual Curves March Food Drive.  This was 68% more than last years total donation!!

 Curves donated to the Yolo County Animal Shelter in January, contributing towels and blankets, detergent, food, toys, beds, and other needed items.

 Curves Davis West Adopts 18 people for Christmas

Curves members adopted eighteen people in the STEAC Program for Christmas and owner Ann Arneson would like to ask readers to join them. STEAC can always use help with cash and food donations. You can contact STEAC by phone at (530) 758 8435, by email at steac@steac.org, or online at http://steac.org. Each year STEAC helps provide food, clothing, and holiday help for over 400 lowincome families, including more than 800 children right here in Davis. Curves put together 45 boxes of food and gifts in November and December 2011 to give back to our community to help those in need.

 Curves member, Marie Risley, created a beautiful fall colored quilt that was raffled to raise money for Thanksgiving food and needs at STEAC, our local food bank.  Over $1200 was donated to help local families have something to be thankful for!

 Curves Davis West was the top fund raising team for the 2011 DAVIS RELAY 4 LIFE.

 Donations

Curves donated $500,000 to Haiti to help with the relief efforts there. Every March they donate millions of pounds of food to help stock local food banks like STEAC. The Curves Davis West Club has donated to these projects as well as the schools, abused & homeless women, STEAC, Foster Families, Progress Ranch, and other local causes.

 Politics/Corporate

There has been some public fallout and controversy, and a loss of members in recent years, when it became public knowledge that the CEO and Founder donated money to some very vocal, antiabortion organizations, although he also donated to a variety of other charity organizations and pregnancy crisis centers in various parts of the country.  A religious convert, he is proud to run his company with the same biblical principles and is outspoken against abortion (http://www.snopes.com/politics/business/curves.asp link).  Curves Davis West lost 1 member because of the previous implication and has grown from 62 to 262 members.  The Curves founders support organizations they personally choose with their personal money, just as all of us do.  That is a freedom we all have.

Please note that Curves is a franchise. A local owner owns the local Curves fitness center. As a franchise the local owner pays to use the Curves name on her fitness center and pays a percentage of the profits to the Curves company.

 History

Davis used to have a Curves in East Davis, but it closed unexpectedly in late June, 2008. A sign on the door said that it closed suddenly due to unforeseen circumstances, and members were encouraged to contact the West Davis or Woodland locations about transferring their membership. However, Woodland also closed in Oct, 2008.  The owner of these two Curves was very ill and has since passed away. He was not able to sell the businesses unfortunately.  Curves Davis West has adopted members from these two clubs as well as Winters Curves. Travelling members from any Curves around the state or world are welcomed to use this Club for free with their travel pass!

Please note that the first comments below do not apply to the current Curves Davis West owner & club. Read further down to get an accurate picture or better yet, just come see for yourself!   The staff are mostly Exercise Bio majors or graduates from UCD and are also well trained by the owner with proven & researched Curves procedures & principles.

Curves Davis West renovated and moved into a bigger space in Westlake Plaza Shopping Center (#106) over Labor Day weekend, 2009.

  See our Gyms page to find out about other gyms in the Davis area.

20041221 14:15:57 nbsp Yes, I sadly quit Curves when I discovered that a good percentage of the profits went toward antiChoice organizations.  Users/BevSykes

20060402 16:49:39 nbsp Being a guy, Ive never been to a Curves before.  I have seen many injuries from patrons at other Curves franchises including shoulder and knee pain that can last for several months.  Ive also known an instructor from Curves who was one of their sucess stories (lost a lot of weight).  Other than that, had no exercise training whatsoever.  Dont keep doing exercises that hurt, no matter WHO recommended them. Users/UsetheForce

20060913 16:46:23 nbsp BORING! This is a good place for women who really are starting from the beginning in terms of fitness (moms, older women mostly). But I found it to be dull and got bored doing the same routine to the same music over and over.  Users/RachelCakes

20070228 01:24:57 nbsp This is the most sexist gym ever.  I thought this would be a great place to go to train my voluptuous body, but they didnt even want me in the gym let alone let me join. Users/BradBenedict

20070430 15:12:50 nbsp Curves is a niche business and its not going to appeal to everyone.  Yes, the company may have some objectionable practices, but I challenge anyone to find a company owned by someone who is completely politically correct for all potential customers, whether its the charities or causes they choose to support or their personal philosophy.  If you like to demonstrate your political and religious beliefs in your commercial activities, Curves may not be for you.  I joined Curves West Davis last summer when it came under new (and the current) ownership after it nearly failed in early 2006, primarily because it was close to my house so I knew I was more likely to go.  The current owner is working very hard to rebuild the club and is very open to suggestions.  Not being Christian, the religious overtones dont motivate me but as a member of a religious minority thats the story of my life.  There is nothing prolife there either. Users/JanCarmikle

20070430 22:49:08 nbsp I only see people doing star steppers and the aerobic blocks, and thats about it. Me and a DDR guy would a try the movements we observed people in curves doing and didnt feel much of a workout, Curves became a running joke because we would see it as we stopped by the new nugget before work Users/StevenDaubert

20070707 22:33:18 nbsp The new owner is very positive.  Many times shes worked out with us.  They are very involved in how you are doing.  You get weighed in once a month and they send postcards.  I love this place.  New owner has replaced some of the music but if you dont like it, take your own.  Great place, I love it!  Im not religious and the owner does not press her religious views on me. I found that to be such a relief.  Users/intangible

20080223 12:29:25 nbsp Does anyone know the monthly fee to join?  Users/HeatherFlood

20080224 16:34:38 nbsp It is $99 every three months. Users/JackkiCox

20080723 15:44:28 nbsp I lost 4 inches in my first 2 months.  The 30 minutes workout is great for busy women, esp moms.  The rate is $44 per month with a 10% discount if you pay for a year.  I dont know where the person thought it was $99 for 3 months.  Also, Curves East Davis just closed so this is the only Curves in Davis. Users/musicamyl

20110404 09:40:32 nbsp Boy, do you folks have it wrong about Curves. Im not talking about the politics, but about the workout and safety of the fitness center. Ive been going over almost two years and have lost (and kept off) 40 pounds, eliminated 2 of the 3 blood pressure medications I was on, and prevented the need for a cholesterol medication. Im stronger, having converted fat to muscle. There are many of us feeling successful and much healthier; thats why women continue. The owner Ann is professional, safety conscious and a strong supporter of the Davis community through food drives, participation in the cancer relay weekend, STEAC, clothing and cosmetic bags for homeless women and much more. I will say the owner of the Curves in east Davis was the exact opposite which is why I left before it closed. Curves Davis West has recently added Zumba to the 30 minute circuit, which has eliminated much of the potential boredom for some. What I like about Curves is that its all women, the machines are hydraulic resistance (you dont have to change each machine as you do with weights), the staff watch and coach you every time, its communityspirited and I can either take my time there and socialize, or get in and out in 30 minutes.

For the month of April women can bring in a $30 bag of groceries for the food drive which waives the $150 membership fee and the first months $44 fee.  You will get a free health evaluation and be taught the circuit routine.  I encourage women of all ages to give Curves in Davis a try!  Users/Laurie

20111226 21:17:05 nbsp I have to agree that while the founder of Curves may be antichoice he has much of his Curves related merchandise (they sell dozens of Curves brand schwag items) manufactured in China and since we all know that countries abortion policy it sort of begs the question. This location is truly great, the owner is top knotch and the rates are affordable. If you are needing to workout and lose weight this is the place because it is not intimidating and will yield results. Woman of all ages with weight to lose will benefit. A few things to keep in mind; you will have someone watching you workout ALL the time (for safety and form) and that can get sometimes be annoying when you know your body and how it feels. Dont do ANY moves that hurt no matter what. The girls at this Curves currently have almost no fitness background apart from Curves, Zumba and CPR training so just take what they say with a grain of salt. The Zumba classes which are totally FUN and are scheduled for a five days a week take precedent over the floor which means if you want to do the circuit you HAVE to do the Zumba class or wait until it ends. you have to donate a dollar to do Zumba but it is SO worth it! Users/MaryCarroll

20120312 11:20:47 nbsp I personally dont care who the CEO/founder supports or doesnt support as long as nobody tries to make me feel like I have to do the same. Since each Curves is an independent franchise, each gym is only as good as the owner/manager and Davis is very lucky to have the best. For those of you who had a negative experience years ago, please come back and give it another try. If youre just watching through the window, it looks very slow and easy.  The fact that you can be in and out in 1/2 hr. might make you think it cant be very effective, but you would be wrong on both counts. Ann and the circuit coaches train every new member on every machine, and everyone is watched to make sure theyre doing everything safely and properly to avoid injuries. The whole atmosphere is very warm, supportive and friendly. The food drives for STEAC, the collections of towels, blankets, food and treats for the animal shelter and the collection of toiletries and personal items for women in shelters for Mothers Day make us all feel like were giving back to the community. And, thanks to Curves, Ive lost 60 lbs. and hope to keep it off by spending 1/2 hr. there every weekday morning. Users/JoyceTrujillo

20120417 09:25:19 nbsp The people here are so friendly and always remember your name. I had to go in there and cancel my membership due to moving, and they were so helpful and made the entire process very easy. The trainers are always there to help you and the owner is always around talking to the customers. I always felt very welcome here. Users/KMcGuire

20120509 15:56:57 nbsp As a person who generally avoided all kinds of physical exercise (though I loved mental exercise like crossword puzzles, cryptograms, etc.), I worried that I doing unseen harm to my body because of my couchpotato ways.  So I thought Id give Curves a try, since it was close to my house and I could go on my own timeschedule.  What a success!  It certainly suits me more than any other form of exercise I have tried, and Ive now stayed with it for about three years, going three times a week.  I appreciate the convenience, the shortness of the workout, the combination of aerobics and weight training (Ive definitely gained strength in my arms), the friendly welcoming atmosphere, and the kindness of the staff.  Theres even information about proper nutrition available, and I learned how to modify my diet with a resulting weight loss of eight pounds.  (I now weigh 120 pounds and am 54 tall.)  The owner even challenges members with weekly quizzes and puzzles, so I get my mental exercise as well!  Barbara Drushell   Users/BarbaraDrushell

20120511 18:07:22 nbsp I found Curves as I sought a place to exercise.  Im sorry to hear the founder of Curves has different values from mine, yet my workouts at Curves Davis West seem quite distant from the corporate headquarters. Sure, I guess ultimately I am contributing to his ability to donate to his political causes which may be in opposition to mine, yet Im not convinced a boycott by me would have much impact on him.  My physical wellbeing as well as the friendships I make at Curves which could open doors to further sharing of my values could ultimately be more effective than my boycott.   The owner, Ann A. has done an outstanding job in creating an open, extremely accepting, warm, friendly environment with a strong emphasis on getting a good solid workout under the supervision of her trainers.  This gym is well utilized daily by women of all sizes, shapes, and I suspect a wide diversity of political and social viewpoints!  We come together to meet our needs to exercise in a safe, clean, inviting environment.  Members enter with a smile and leave with a smile and a wave! All clients are welcomed, introduced, and acknowledged daily by first names while they exercise, building a wonderful collegiality that keeps you coming back.  (Ann also stays on top of her members so that postcards or phone calls come in if you havent been in for a while!) Ann also puts a high value on giving back to her community through the various service projects she organizes and encourages her members to participate in, STEAC, Womens Shelter, womens health awareness, etc. If youre looking for a friendly, supportive place for women to exercise together, check out Curves Davis West! SCole
 Users/charlottecole

20120709 15:02:20 nbsp yet my workouts at Curves Davis West seem quite distant from the corporate headquarters.   What an odd statement. You support Curves, you support Curvesend of. Curves is a business with countless items made in China and I will not support CEO Gary Heavein and Chinas liquidation of the worlds last Rhinos, Tigers (and God knows what else), thats not Christianity thats corporate greed under the guise of just locals helping locals. Truly happy switched Peak Performance. The girls who work at Curves are appallingly rude and whats the point of being jazzed about fitness with someone staring at you constantly? Great that the bulk of their charity work goes to help folks in Europe though.... Users/MaryCarroll

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      module phys_constants

      real(kind=8), parameter :: h=6.626d-27,k=1.38d-16, &
       c=2.99792458d10,e=4.8032d-10,G=6.67d-8, &
       m=9.10938188d-28,mp=1.67262158d-24,pi=acos(-1d0),c2=c*c,sigb=5.6704d-5, &
       msun=1.998d33, sigt=6.6523d-25

      contains

        subroutine phys_constants_dummy()
        end subroutine phys_constants_dummy

      end module

max_stars_repo_name: HerculesJack/grtrans

max_issues_repo_name: RAnantua/grtrans

max_forks_repo_name: RAnantua/grtrans


      program main
!----------------------------------------------------------------
!     trajectory guided gaussian basis
!----------------------------------------------------------------------
      use cdat

      implicit real*8(a-h,o-z)

!      real*8, dimension(nmax) :: ke,v,x,u,p,w,du,qp,dv,ddv,s,q,s0

      real*8,allocatable :: ke(:),v(:),x(:),u(:),p(:),w(:)
      real*8,allocatable :: du(:),qp(:),dv(:),ddv(:),s(:)

      complex*16,dimension(:,:),allocatable :: h,mat
      complex*16, allocatable :: c(:),c0(:),dc(:),c1(:)
      complex (kind = 8), allocatable ::  cold(:), cnew(:)

!      complex*16, dimension(nmax,nmax) :: mat

!      complex*16, dimension(nmax) :: c,c0,dc

     	real*8 :: f(3),ki,ncm

      real*8 :: gasdev

      complex*16 :: z0,z1,z2,zn,z,cor

!      common/wave/q0,p0,a0,sigma
!      common/grid/xmin,xmax,dx,np

      open(100,file='energy.dat')
      open(101,file='phase.dat')
      open(102,file='aver')
      open(103,file='wft.dat')
      open(104,file='wf0.dat')
      open(106,file='q.dat')
      open(105,file='p.dat')
      open(107,file='c.dat')
      open(108,file='basis.dat')
!      open(109,file='basisn.dat')
      open(110,file='cf0.dat')
      open(111,file='norm.dat')
      open(112,file='force')
      open(113,file='xc')
      open(114,file='temp.dat',status='unknown',action='write')
      open(115,file='cor.dat', status='unknown',action='write')
      open(116,file='qforce.dat', status='unknown',action='write')
      open(117,file = 'pes.dat')

      call init_pars()

      ep = 1d-2
      Nb = Ntraj
      Np = 800
      af = 0.1

      write(*,1002) ntraj
1002  format(/'1D QTM code with Gassian Basis '//,  &
              'analytic expression for matrix elements'//, &
              'system: Morse',//        &
              'number of basis:', i4//)

      allocate(ke(Ntraj),v(Ntraj),x(Ntraj),p(Ntraj),s(Ntraj), &
               w(Ntraj),du(Ntraj),u(Ntraj),dv(Ntraj),ddv(Ntraj))

      allocate(mat(Ntraj,Ntraj))
      allocate(c(Ntraj),c0(Ntraj),h(Ntraj,Ntraj),dc(Ntraj),c1(ntraj))

      allocate(cold(ntraj),cnew(ntraj))

      pow = 6d0

!     grid size over whole range

      xmin = -10d0
      xmax = 10d0
      dx = (xmax-xmin)/dble(np-1)

! --- plot the potential
      do i=1,np
        q = xmin + dx*(i-1)
        call pot(q,vpot,addv)
        write(117,1000) q,vpot
      enddo

! --- grid to sample initial points
      xl = q0-1.d0*dsqrt(pow/a0)
      xr = q0+1.d0*dsqrt(pow/a0)
      gx = (xr-xl)/dble(nb-1)

      do i=1,nb
!          x(i) = gasdev(idum1)
!          x(i) = x(i)/dsqrt(2d0*a0)+q0
          p(i) = 0d0
          x(i) = xl+(i-1)*gx
          w(i) = exp(-2d0*a0*(x(i)-q0)**2)*dsqrt(2d0*a0/pi)*gx
      enddo

! --- action function

      s = 0d0

! --- define gaussian width as half of bin size
!      al = 1d0/2d0/sigma**2
      sigma = sqrt(1d0/2d0/al)

!      xint = 1.8d0*sigma

!     test distribution
      an = 0d0
      a1 = 0d0
      do i=1,Ntraj
        an = an+w(i)
        a1 = a1+w(i)*x(i)
      enddo
      write(*,*) an,a1
      write(*,*) xl,xr,gx

!     print out the initial conditions
      write(*,1001) sigma,a0,Ntraj,kmax,dt,kout,am,p0,q0
1001  format('Initial Conditions'//, &
            'width = ', f10.6/, &
            'a0    = ', f10.6/, &
            'Ntraj = ', i6/ ,   &
            'Kmax  = ', i6/ ,   &
            'dt    = ', f10.6/, &
            'kout  = ', i6/ ,   &
            'Mass  = ', f10.6/ ,&
            'p0    = ', f10.6/ ,&
            'q0    = ', f10.6//)

!--------initial coeff---------------------
      write(*,1011)
1011  format(/'Computing inital coefficients'/)

      call proj(nb,x,p,s,c0)

      c = c0
      do i=1,nb
        write(110,1000) x(i),abs(c(i)),c(i)
      enddo
      write(6,1012) x
1012  format('initial sampling position = '//,20(f5.2,1x)/)

      write(6,1014) p
1014  format('initial sampling momentum = '//,20(f5.2,1x)/)

      write(6,1015) c
1015  format('initial expansion coeff   = '//,20('(',f5.2,',',f5.2,')',1x)/)

!------check normalization--------------
      z1 = (0d0,0d0)
      do i=1,np
        xi = xmin+dx*(i-1)
        call psi(xi,nb,c,x,s,z0)
        z1 = z1+conjg(z0)*z0*dx
        if(abs(z0).gt.1d-3) write(104,1000) xi,abs(z0)**2
      enddo

      call norm(nb,x,p,c,s,anrm)
      write(*,1013)  anrm,real(z1)
1013  format('initial normalization =', 2f12.5/)

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

!     print out basis functions
      do i=1,Np
        xi = xmin+dx*dble(i-1)
        z0 = dsqrt(dsqrt(al/pi))*exp(-al/2d0*(xi-x(1))**2+  &
             im*p(1)*(xi-x(1)))
        zn = dsqrt(dsqrt(al/pi))*exp(-al/2d0*(xi-x(2))**2+ &
             im*p(nb)*(xi-x(2)))

       if(abs(z0)>1d-3 .or. abs(zn)>1d-3) write(108,1000) xi,abs(z0)**2,abs(zn)**2

!        if(abs(z0).gt.1d-3) write(108,1000) xi,abs(z0)**2
!        if(abs(zn).gt.1d-3) write(109,1000) xi,abs(zn)**2
      enddo

! --- initial force field

!      call derivs(x,Ntraj,v,dv,ddv)
!      call qpot(order,nb,ddv,c,am,p,x,s,w,u,du)

!      do i=1,nb
!        write(112,1000) x(i),v(i),dv(i),ddv(i),du(i)
!      enddo

!      call ham(am,w,c,nb,x,p,s,u,du,v,ddv,h)
!     call increm(nb,mat,h,c,dc)

!-----output at t=0-------------------
      t = 0.0
      write(107,1000) t,(abs(c(i)),i=1,Ntraj)
      dt = dt/kout
      dt2 = dt/2d0
      af = 0.0 !friction constant

      call ham0(c,nb,x,p,s,dc)

!      c1 = c
      cold = c
      c = c + dc*dt
!      x = x + p*dt/am

!-----begin the time step loop----------
      time: do kt = 1,kmax

!        print *,'time step', kt

        do 13 kk=1,kout

          t = t+dt

!          call norm(nb,x,p,c,mat,an)

!          c = c + dc*dt2

          call mom(nb,x,p,c,s)

          x = x + p*dt/am

!-------------update dc/dt--------------------
!          call derivs(nb,x,v,dv,ddv)

        call ham0(c,nb,x,p,s,dc)

        cnew = cold + dc*2d0*dt
        cold = c
        c = cnew

! --------- friction
          call norm(nb,x,p,c,s,anm)
          c = c/dsqrt(anm)

13      enddo

!        if(k == kmax) then
!        do i=1,nb
!          write(116,1000) x(i),du(i)
!        enddo
!        endif

! ----- correlation function

        call corr(nb,x,p,c,s,cor)

        write(115,1000) 2d0*t,cor,abs(cor)

        write(107,1000) t,(abs(c(i)),i=1,nb)

!       print out trajectories in phase space
        write(101,1000) (x(i),i=1,nb),(p(i),i=1,nb)

        write(105,1000) t,(p(i),i=1,Nb)
        write(106,1000) t,(x(i),i=1,nb)
!       calculate the total energy
        enk = 0d0
        env = 0d0
        enq = 0d0

!        do i=1,Nb
!          env = env+v(i)*w(i)
!          enk = enk+p(i)**2/2d0/am*w(i)
!          enq = enq+u(i)*w(i)
!        enddo

!        en = enk+env+enq

!        write(100,1000) t,enk,env,enq,en
!        call flush(100)

!-------------average position --------------
!        call aver(Ntraj,x,p,c,xav)
!        write(102,1000) t,xav

      end do time

!--------print x,c--------------
        do i=1,nb
          write(113,1000) x(i),abs(c(i))
        enddo
!---------check matrix singularity--------------
!        do i=1,nb
!          write(*,1018) (mat(i,j),j=1,nb)
!        enddo
!        call eigen(nb,mat)

!1018    format(20(f6.4,1x))

!     print out final wavefunction
      prob = 0d0
      do i=1,Np
        xi = xmin + dx*(i-1)
        call psi(xi,nb,c,x,s,z0)
        if(abs(z0).gt.1d-3) write(103,1000) xi,abs(z0)**2
        prob = prob+abs(z0)**2*dx
      enddo

      write(*,*) 'Final prob = ',prob
      write(114,*) t,gx
      write(114,*) kmax,',',dt,',',kout
      write(114,*) am
      write(114,*) idum1
      write(114,*) order
      write(114,*) nb
      write(114,*) al,beta

      write(114,*) (x(i),i=1,nb)
!      write(114,*) (p(i),i=1,nb)
!      write(114,*) (s(i),i=1,nb)
      write(114,*) (c(i),i=1,nb)

      close(111)
      close(112)
      close(113)
      close(114)
      close(115)

1000 format(2000(e14.7,1x))
      end program main

! ---- update momentum

      subroutine mom(nb,x,p,c,s)

      use cdat, only : al,pi,beta

      implicit real*8(a-h,o-z)

      real*8 x(nb),p(nb),s(nb)
      complex*16 c(nb),z,z0,dz

! ----- Gaussian convolution

!      beta = 16d0

      an0 = dsqrt(dsqrt(al/pi))
      an0 = an0*dsqrt(beta/(beta+al))
      alfa = beta*al/(beta+al)

      do j=1,nb
        z = (0d0,0d0)
        dz = (0d0,0d0)
        q = x(j)
        do i=1,nb
          y = q-x(i)
          z0 = c(i)*exp(-0.5*alfa*y**2)*an0
          z = z + z0
          dz = dz - alfa*y*z0
        enddo
        p(j) = imag(dz/z)
      enddo

      return
      end subroutine
!------------------------------------------------

! --- autocorrelation
!      double complex function corr(nb,c,q,p,s)
!
!      use cdat
!
!      implicit real*8 (a-h,o-z)
!
!      real*8, intent(in),dimension(nb) :: q,p,s
!      complex*16, intent(in) :: c(nb)
!
!      complex*16 :: mat(nb,nb),z,z0,c0(nb)
!!      common/wave/q0,p0,a0,sigma
!
!      a = al
!      alfa = 2d0*a0
!      as=a+alfa
!      av=a*alfa/as
!      an=dsqrt(2d0*dsqrt(av/as))
!      an0=dsqrt(dsqrt(a/pi))
!      an2=dsqrt(dsqrt(alfa/pi))
!!---- compute c0_k = <psi0|g_k>
!      do j=1,nb
!        dq=q(j)-q0
!!        dp=p(j)-p0
!        dp = -p0
!        d2=-0.5d0*av*dq*dq
!        d0=-0.5d0/as*dp*dp
!        d1=(a*p0)/as*dq
!        c0(j)=exp(d2+d0+im*d1)*an
!      enddo
!
!      corr = (0d0,0d0)
!      do i=1,nb
!        corr = corr+conjg(c0(i))*exp(im*s(i))*c(i)
!      enddo
!
!      end function

! --- compute wavefunction value at point x

      subroutine psi(xx,nb,c,q,s,z)

      use cdat, only : im,pi,al

      implicit real*8(a-h,o-z)
      real*8 ,intent(in) :: xx
      integer*4 :: nb
      real*8, intent(in), dimension(nb) :: q,s
      complex*16,intent(in) :: c(nb)

      complex*16, intent(out) :: z

      anrm = dsqrt(dsqrt(al/pi))

      z = (0d0, 0d0)
      do j=1,nb
        z = z + c(j)*exp(-al/2d0*(xx-q(j))**2+im*s(j))*anrm
      enddo

      end subroutine

! --- compute non-classical momemtum

!      double precision function ncm(x,nb,c,q,p,s)

!      use cdat, only : im,pi,al

!      implicit real*8(a-h,o-z)
!      real*8 ,intent(in) :: x
!      integer*4 :: nb
!      real*8, intent(in),dimension(nb) :: q,p,s
!      complex*16 :: c(nb),psi,dpsi,z0,z1,gi,gj

!      common/wave/q0,p0,a0,sigma
!      anrm = dsqrt(dsqrt(al/pi))

!      ncm = 0d0
!      do i=1,nb
!        ncm = ncm+abs(c(i))**2*anrm**2*exp(-al*(x-q(i))**2)
!      enddo

!      h = 1d-4
!      z0 = psi(x,nb,c,q,p,s)
!      xx = x+h
!      z1 = psi(xx,nb,c,q,p,s)
!      xx = x-h
!      z2 = psi(xx,nb,c,q,p,s)
!      ncm = (abs(z2)-abs(z1))/(2d0*h*abs(z0))

!      den = conjg(z0)*z0
!      den1 = conjg(z1)*z1
!      ncm = (den1-den0)/(2d0*h*den0)

!      do i=2,nb
!        do j=1,i-1
!          dq = q(i)-q(j)
!          dp = p(i)-p(j)
!          xs = (q(i)+q(j))/2d0
!          ps = (p(i)+p(j))/2d0
!          ncm = ncm+anrm**2*real(conjg(c(i))*c(j)*(-al*(2d0*x-q(i)-q(j))-im*dp) &
!                *exp(-al*(x-xs)-im*dp*(x-xs)+im*ps*dq+im*(s(j)-s(i))))
!        enddo
!      enddo
!      ncm = ncm/den
!      do i=1,nb
!        do j=1,nb
!          dp = p(i)-p(j)
!          gi = anrm*exp(-al/2d0*(x-q(i))**2+im*p(i)*(x-q(i))+im*s(i))
!          gj = anrm*exp(-al/2d0*(x-q(j))**2+im*p(j)*(x-q(j))+im*s(j))
!          ncm = ncm+conjg(c(i))*c(j)*(-al*(2d0*x-q(i)-q(j))-im*dp)*conjg(gi)*gj
!        enddo
!      enddo

!      ncm = ncm/(2d0*den)

! --- analytical way to compute r(x)
!      z0 = psi(x,nb,c,q,p,s)
!      z1 = dpsi(x,nb,c,q,p,s)

!      ncm = real(conjg(z0)*z1/abs(z0)**2)

!      end function

      subroutine dpsi(x,nb,c,q,s,z)

      use cdat, only : im,pi,al

      implicit real*8(a-h,o-z)

      real*8 ,intent(in) :: x
      integer*4 :: nb
      real*8, intent(in), dimension(nb) :: q,s
      complex*16,intent(in) :: c(nb)
      complex*16 :: gj

      complex*16, intent(out) :: z

      anrm = dsqrt(dsqrt(al/pi))

      z = (0d0, 0d0)
      do j=1,nb
        gj = exp(-al/2d0*(x-q(j))**2+im*s(j))*anrm
        z = z + c(j)*(-al*(x-q(j)))*gj
      enddo

      end subroutine

!------normalizati-----------------------------
      subroutine norm(nb,x,p,c,s,anrm)

      use cdat, only : pi, al, im

      implicit real*8(a-h,o-z)

      integer*4, intent(in) :: nb

      real*8, dimension(nb) :: x,p,s

      complex*16, intent(in) :: c(nb)

      complex*16 :: mat(nb,nb),z

      real*8, intent(out) :: anrm

      z = (0d0,0d0)

      call multi(nb,x,p,s,mat)

      do j=1,nb
        do k=1,nb
          z = z + conjg(c(j))*mat(j,k)*c(k)
        enddo
      enddo

      anrm = abs(z)

      return
      end subroutine

      subroutine corr(nb,x,p,c,s,cor)

      use cdat, only : pi, al, im

      implicit real*8(a-h,o-z)

      integer*4, intent(in) :: nb

      real*8, dimension(nb) :: x,p,s

      complex*16, intent(in) :: c(nb)

      complex*16 :: mat(nb,nb),z,cor

      z = (0d0,0d0)

      call multi(nb,x,p,s,mat)

      do j=1,nb
        do k=1,nb
          z = z + c(j)*mat(j,k)*c(k)
        enddo
      enddo

      cor = z

      return
      end subroutine

!     ----------------------------------------------
!     average
!     ---------------------------------------------
!      subroutine aver(q,p,s,c,xav)
!      use cdat
!
!      implicit real*8(a-h,o-z)
!      integer*4,intent(in) :: Ntraj
!      real*8,intent(in)::q(Ntraj),p(Ntraj),s(ntraj)
!      complex*16, intent(in) :: c(Ntraj)
!      real*8,intent(out) :: xav
!      complex*16 :: psi,z
!!      common/wave/q0,p0,a0,sigma
!      common/grid/xmin,xmax,dx,np
!
!      z = (0d0,0d0)
!      do i=1,Np
!        xi = xmin+dx*(i-1)
!        psi = (0d0,0d0)
!        do j=1,ntraj
!          psi = psi+c(j)*exp(-al/2d0*(xi-q(j))**2)*dsqrt(dsqrt(al/pi))
!        enddo
!        z = z+conjg(psi)*psi*dx*xi
!      enddo
!
!      xav = real(z)
!
!      return
!      end subroutine
!     -----------------------------------
!     gassian wavepackets overlap matrix, analytic solution
!     -----------------------------------
      subroutine multi(nb,q,p,s,mat)

      use cdat, only : im,pi,al

      implicit real*8(a-h,o-z)
      integer*4, intent(in) :: nb
      real*8   , intent(in) :: q(nb),p(nb),s(nb)
      complex*16,intent(out) :: mat(nb,nb)
      complex*16 :: z,z0
!      common/wave/q0,p0,a0,sigma

      anrm = dsqrt(dsqrt(al/pi))
      do j=1,nb
        do k=1,j
          dq = q(j)-q(k)
          z = dq/2d0
          z0 = exp(-al*z*conjg(z))
          mat(j,k) = z0*exp(im*(s(k)-s(j)))
        enddo
      enddo

      do j=2,nb
        do k=1,j-1
          mat(k,j) = conjg(mat(j,k))
        enddo
      enddo

      return
      end subroutine
! --------------------------------------------
      subroutine pot(x,v,ddv)
      implicit real*8(a-h,o-z)

      ipot = 2

      if (ipot == 0) then
        ep = 0.0
        v = x**2/2d0 + ep*x**4/4.0
        ddv = 1d0 + 3.0*x**2*ep

      elseif(ipot == 1) then

        eta = 1.3544d0
        v = 1d0/16d0/eta*x**4-0.5d0*x**2
        ddv = 3d0/4d0/eta*x**2-1d0
      elseif(ipot == 2) then

        D = 16d0
!        a = 1.3624d0
        a = 1d0
        v = D/cosh(a*x)**2
        ddv = 6d0*D*a**2*sinh(a*x)**2/cosh(a*x)**4-2d0*D*a**2/cosh(a*x)**2
      endif

      return
      end subroutine

      subroutine ham0(c,nb,q,p,s,dc)

      use cdat

      implicit real*8(a-h,o-z)

      integer*4, intent(in) :: nb
      complex*16, dimension(nb), intent(in) :: c
      real*8, intent(in), dimension(nb) :: q,p,s

      complex*16,intent(out) :: dc(nb)

      real*8 :: v,dv,ddv

      complex*16, dimension(nb,nb) :: h,mat
      real*8 :: aux(nb,nb)

      complex*16 :: z,z0,d0,d1,d2,b(nb)
      real*8 b1(nb,2)

      call multi(nb,q,p,s,mat)

      anrm = dsqrt(dsqrt(al/pi))
      eta = 1.3544d0

      do j=1,nb
        do k=1,nb
          dq = q(k)-q(j)
          z = dq/2d0
          z0 = exp(-al*dq**2/4d0)
          d0 = al/(4d0*am)*(1d0-al/2d0*dq**2)
          y = (q(j)+q(k))/2d0

! --- local harmonic approximation
          call pot(y,v,ddv)
          d1 = v + ddv/(4d0*al)

! --- analytical potential integration for double-well
!          d1 = 1d0/(16d0*eta)*(y**4 + 3d0*y**2/al + 3d0/(4d0*al**2)) - 0.5d0*(y**2+1d0/(2d0*al))

          d2 = im*p(k)*al/(2d0*am)*dq
          h(j,k) = z0*(d0+d1+d2)
        enddo
      enddo

      b = matmul(h,c)

      do i=1,nb
        b1(i,1) = real(b(i))
        b1(i,2) = imag(b(i))
      enddo

!      call inverse(mat,mativs,nb)

!-------save the overlap matrix-----------
      do i=1,nb
      do j=1,nb
        aux(i,j) = real(mat(i,j))
      enddo
      enddo

      call dposv('U', nb, 2, aux, nb, b1, nb, INFO)
      if(INFO .ne. 0) then
        write(*,*) 'dC/dt: Matrix fails, INFO =',INFO
        stop
      endif
!------compute increment dc------------
      do j=1,nb
        dc(j) = -im*(b1(j,1)+im*b1(j,2))
      enddo

      return
      end subroutine

!     -------------------------------------------
!     hamiltonian matrix (H-id/dt) with gaussian
!     basis
!     -------------------------------------------
      subroutine ham(c0,c,nb,q,p,s)

      use cdat, only:pi,im,al,dt

      implicit real*8(a-h,o-z)

      integer*4, intent(in) :: nb
      complex*16, dimension(nb) :: c,c0
      real*8, intent(in), dimension(nb) :: q,p,s

!      complex*16,intent(out) :: dc(nb)

      real*8 :: v,dv,ddv

      complex*16, dimension(nb,nb) :: h,mat
      real*8 :: aux(nb,nb)

      complex*16 :: z,z0,d0,d1,d2,b(nb)
      real*8 b1(nb,2)

!      common/wave/q0,p0,a0,sigma

      call multi(nb,q,p,s,mat)

      anrm = dsqrt(dsqrt(al/pi))

      do j=1,nb
        do k=1,nb
          dq = q(k)-q(j)
          z = dq/2d0
          z0 = exp(-al*dq**2/4d0)
!          d0 = z0*(v(k)-(p(k)**2-al)/2d0/am) ! without action
          d0 = al/(4d0*am)*(1d0-al/2d0*dq**2)
!          d1 = v(k) + dv(k)*z + 0.5*ddv(k)*(z**2 + 1d0/(2d0*al))
          y = (q(j)+q(k))/2d0
          call pot(y,v,ddv)
          d1 = v + ddv/(4d0*al)
          d2 = im*p(k)*al/(2d0*am)*dq
          h(j,k) = z0*(d0+d1+d2)

        enddo
      enddo

!      print *,'h = ', h(1,2),h(3,4)
!      print *,'mat =', mat(1,2),mat(3,4)

!      b = (0d0,0d0)
!      do i=1,nb
!        do j=1,nb
!          b(i) = b(i) + h(i,j)*c(j)
!        enddo
!      enddo
      b = -2d0*dt*im*matmul(h,c) + matmul(mat,c0)

      do i=1,nb
        b1(i,1) = real(b(i))
        b1(i,2) = imag(b(i))
      enddo

!      call inverse(mat,mativs,nb)

!-------save the overlap matrix-----------
      do i=1,nb
      do j=1,nb
        aux(i,j) = real(mat(i,j))
      enddo
      enddo

      call dposv('U', nb, 2, aux, nb, b1, nb, INFO)
      if(INFO .ne. 0) then
        write(*,*) 'dC/dt: Matrix fails, INFO =',INFO
        stop
      endif
!------compute increment dc------------

      c0 = c ! update c0

      do j=1,nb
        c(j) = -im*(b1(j,1)+im*b1(j,2))
      enddo

      return
      end subroutine

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

      subroutine increm(nb,mat,h,c,dc)
      use cdat, only:im,pi
      implicit real*8(a-h,o-z)
      complex*16,intent(in) :: mat(nb,nb),h(nb,nb),c(nb)
      complex*16,intent(out) :: dc(nb)
      complex*16 b(nb),aux(nb,nb)
!---------H*c---------------
!      do i=1,Nb
!        b(i) = (0d0,0d0)
!        do j=1,Nb
!          b(i) = b(i)+h(i,j)*c(j)
!        enddo
!      enddo
      b = matmul(h,c)
!      b = -im*b

!-------save the overlap matrix-----------
      aux = mat
!      do i=1,Nb
!        write(*,1023) (aux(i,j),j=1,nb)
!      enddo
!      write(*,*) (b(i),i=1,Nb)
!1023  format(100(f6.4,1x))
      call zposv('U', Nb, 1, aux, nb, b, nb, INFO)
      if(INFO .ne. 0) then
        write(*,*) 'dC/dt: Matrix fails, INFO =',INFO
        stop
      endif
!------compute increment dc------------
      do j=1,nb
        dc(j) = -im*b(j)
      enddo
!      dc = b

      return
      end subroutine
!     ----------------------------------------------
!     initial coefficients before gaussian basis
!     ----------------------------------------------
      subroutine proj(nb,x,p,s,c)

      use cdat

      implicit real*8(a-h,o-z)
      integer*4,intent(in) :: nb
      real*8,intent(in),dimension(nb) :: x,p,s

      complex*16,intent(out) :: c(nb)

      complex*16 :: c0(nb)
      complex*16 :: z,z0,mat(nb,nb),psi0,aux(nb,nb)

!      common/wave/q0,p0,a0,sigma
      a = al
      alfa = 2d0*a0
      as=a+alfa
      av=a*alfa/as
      an=dsqrt(2d0*dsqrt(av/as))
      an0=dsqrt(dsqrt(a/pi))
      an2=dsqrt(dsqrt(alfa/pi))
!------compute c_k=<g_k|psi0>---------
      do j=1,nb
        dq=x(j)-q0
        dp=p(j)-p0
        d2=-0.5d0*av*dq*dq
        d0=-0.5d0/as*dp*dp
        d1=(alfa*p(j)+a*p0)/as*dq
        c0(j)=exp(d2+d0+im*d1)*an
      enddo

!      write(*,*) 'overlap',mat(1,1),mat(1,2),mat(1,3)
!     intial wavepacket is also a gaussian (p0,q0) with the same width
!     as basis
!      do j=1,nb
!        dq = q(j)-q0
!        z = dq/2d0-im/2d0/al*(p(j)-p0)
!        z0 = exp(-al*z*conjg(z)+im/2d0*(p(j)+p0)*dq)
!        c0(j) = z0*dsqrt(dsqrt(al/pi))
!      enddo

!     a0 not equal to al

!      do j=1,nb
!        z = (0d0,0d0)
!        do i=1,Np
!          x = xmin+dx*(i-1)
!          z0 = exp(-al/2d0*(x-q(j))**2+im*p(j)*(x-q(j)))
!          psi0 = exp(-a0/2d0*(x-q0)**2+im*p0*(x-q0))
!          z = z+conjg(z0)*psi0*dx*dsqrt(dsqrt(a0*al)/pi)
!        enddo
!        c0(j) = z
!      enddo

!------save overlap matrix------------------

      call multi(nb,x,p,s,mat)

      aux = mat

!     solve matrix equation, Mc = b
!      do i=1,nb
!        write(*,1009) (mat(i,j),j=1,Nb)
!      enddo
!1009  format(20('(',f5.2,',',1x,f5.2,')',2x))

      write(6,1010) c0
1010  format('c = ',20('(',f5.2,',',1x,f5.2,')',2x))

      call zposv('U', Nb, 1, aux, nb, c0, nb, INFO)

      IF( INFO.ne.0 ) THEN
         WRITE(*,1222) INFO
1222     format('INFO = ', i6/ ,'Initial projection fails.'/)
         STOP
      END IF

      do i=1,nb
        c(i) = c0(i)
      enddo

      return
      end subroutine

!      subroutine reproj(nb,x,p,s,c)
!      use cdat
!
!      implicit real*8(a-h,o-z)
!      integer*4,intent(in) :: nb
!      real*8,intent(in),dimension(nmax) :: x,p,s
!
!      complex*16,intent(inout) :: c(nmax)
!
!      complex*16 :: c0(nb)
!
!      complex*16 :: z,z0,mat(Nb,Nb),psi0,aux(nb,nb)
!
!      do i=1,nb
!        c0(i) = c(i)
!      enddo
!
!      call multi(nb,x,p,s,mat)
!
!      aux = mat
!
!      call zposv('U', Nb, 1, aux, nb, c0, nb, INFO)
!
!      IF( INFO.ne.0 ) THEN
!         WRITE(*,1222) INFO
!1222     format('INFO = ', i6/ ,'projection fails.'/)
!         STOP
!      END IF
!
!      do i=1,nb
!        c(i) = c0(i)
!      enddo
!
!      return
!      end subroutine


text:
FUNCTION random_Poisson(mu, first) RESULT(ran_Poisson)
!**********************************************************************

!     Translated to Fortran 90 by Alan Miller from:
!                                    RANLIB
!
!           Library of Fortran Routines for Random Number Generation
!
!                           Compiled and Written by:
!
!                                Barry W. Brown
!                                 James Lovato
!
!                    Department of Biomathematics, Box 237
!                    The University of Texas, M.D. Anderson Cancer Center
!                    1515 Holcombe Boulevard
!                    Houston, TX      77030
!
! This work was supported by grant CA-16672 from the National Cancer Institute.

!                 GENerate POIsson random deviate Function

!  Generates a single random deviate from a Poisson
!  distribution with mean mu.

!                           Arguments

!  mu --> The mean of the Poisson distribution from which
!         a random deviate is to be generated.
!                           REAL mu

!                           Method

!  For details see:

!      Ahrens, J.H. and Dieter, U.
!      Computer Generation of Poisson Deviates
!      From Modified Normal Distributions.
!      ACM Trans. Math. Software, 8, 2 (June 1982), 163-179

!  MUPREV = PREVIOUS MU, MUOLD = MU AT LAST EXECUTION OF STEP P OR B.
!  TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
!  COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K) - DEL

!  SEPARATION OF CASES A AND B

IMPLICIT NONE

!     .. Scalar Arguments ..
REAL, INTENT(IN)    :: mu
LOGICAL, INTENT(IN) :: first
INTEGER             :: ran_Poisson

INTERFACE
  FUNCTION random_normal() RESULT (ran_norm)
    IMPLICIT NONE
    REAL :: ran_norm
  END FUNCTION random_normal
END INTERFACE

INTERFACE
  FUNCTION random_exponential() RESULT (ran_exp)
    IMPLICIT NONE
    REAL :: ran_exp
  END FUNCTION random_exponential
END INTERFACE

!     .. Local Scalars ..
REAL    :: a0 = -.5, a1 = .3333333, a2 = -.2500068, a3 = .2000118,        &
           a4 = -.1661269, a5 = .1421878, a6 = -.1384794, a7 = .1250060,  &
           b1, b2, c, c0, c1, c2, c3, d, del, difmuk, e, fk, fx, fy, g,   &
           omega, p, p0, px, py, q, s, t, u, v, x, xx
INTEGER :: j, k, kflag, l, m
!     ..
!     .. Local Arrays ..
REAL    :: fact(10) = (/ 1., 1., 2., 6., 24., 120., 720., 5040., 40320., 362880. /)
REAL    :: pp(35)
!     ..
!     ..
!     .. Executable Statements ..
IF (.NOT. first) GO TO 10
IF (mu < 10.0) GO TO 120

!     C A S E  A. (RECALCULATION OF S, D, L IF MU HAS CHANGED)

s = SQRT(mu)
d = 6.0*mu*mu

!     THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL PROBABILITIES FK
!     WHENEVER K >= M(MU).
!     L = INT(MU-1.1484) IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .

l = INT(mu-1.1484)

!     STEP N. NORMAL SAMPLE - random_normal() FOR STANDARD NORMAL DEVIATE

10 g = mu + s*random_normal()
IF (g < 0.0) GO TO 20
ran_Poisson = INT(g)

!     STEP I. IMMEDIATE ACCEPTANCE IF ran_Poisson IS LARGE ENOUGH

IF (ran_Poisson.GE.l) RETURN

!     STEP S. SQUEEZE ACCEPTANCE - SAMPLE U

fk = REAL(ran_Poisson)
difmuk = mu - fk
CALL RANDOM_NUMBER(u)
IF (d*u >= difmuk*difmuk*difmuk) RETURN

!     STEP P. PREPARATIONS FOR STEPS Q AND H.
!             (RECALCULATIONS OF PARAMETERS IF NECESSARY)
!             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
!             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
!             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
!             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.

20 IF (.NOT. first) GO TO 30
omega = .3989423/s
b1 = .4166667E-1 / mu
b2 = .3*b1*b1
c3 = .1428571*b1*b2
c2 = b2 - 15.*c3
c1 = b1 - 6.*b2 + 45.*c3
c0 = 1. - b1 + 3.*b2 - 15.*c3
c = .1069/mu
30 IF (g < 0.0) GO TO 50

!             'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)

kflag = 0
GO TO 70

!     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)

40 IF (fy-u*fy <= py*EXP(px-fx)) RETURN

!     STEP E. EXPONENTIAL SAMPLE - random_expon() FOR STANDARD EXPONENTIAL
!             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
!             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)

50 e = random_exponential()
CALL RANDOM_NUMBER(u)
u = u + u - 1.0
t = 1.8 + SIGN(e, u)
IF (t <= (-.6744)) GO TO 50
ran_Poisson = INT(mu+s*t)
fk = REAL(ran_Poisson)
difmuk = mu - fk

!             'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)

kflag = 1
GO TO 70

!     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)

60 IF (c*ABS(u) > py*EXP(px+e) - fy*EXP(fx+e)) GO TO 50
RETURN

!     STEP F. 'SUBROUTINE' F. CALCULATION OF PX, PY, FX, FY.
!             CASE ran_Poisson < 10 USES FACTORIALS FROM TABLE FACT

70 IF (ran_Poisson >= 10) GO TO 80
px = -mu
py = mu**ran_Poisson / fact(ran_Poisson + 1)
GO TO 110

!      CASE ran_Poisson >= 10 USES POLYNOMIAL APPROXIMATION
!      A0-A7 FOR ACCURACY WHEN ADVISABLE
!      .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)

80 del = .8333333E-1 / fk
del = del - 4.8*del*del*del
v = difmuk/fk
IF (ABS(v) <= 0.25) GO TO 90
px = fk*LOG(1.0+v) - difmuk - del
GO TO 100

90 px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) - del
100 py = .3989423/SQRT(fk)
110 x = (0.5-difmuk)/s
xx = x*x
fx = -0.5*xx
fy = omega* (((c3*xx+c2)*xx+c1)*xx+c0)
IF (kflag > 0) THEN
  GO TO 60
ELSE
  GO TO 40
END IF

!     C A S E  B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)

120 IF (first) THEN
  m = MAX(1, INT(mu))
  l = 0
  p = EXP(-mu)
  q = p
  p0 = p
END IF

!     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD

130 CALL RANDOM_NUMBER(u)
ran_Poisson = 0
IF (u <= p0) RETURN

!     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
!             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
!             (0.458 = PP(9) FOR MU = 10)

IF (l == 0) GO TO 150
j = 1
IF (u > 0.458) j = MIN(l, m)
DO k = j, l
  IF (u <= pp(k)) GO TO 180
END DO
IF (l == 35) GO TO 130

!     STEP C. CREATION OF NEW POISSON PROBABILITIES P
!             AND THEIR CUMULATIVES Q = PP(K)

150 l = l + 1
DO k = l, 35
  p = p*mu/REAL(k)
  q = q + p
  pp(k) = q
  IF (u <= q) GO TO 170
END DO
l = 35
GO TO 130

170 l = k
180 ran_Poisson = k
RETURN

END FUNCTION random_Poisson

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
! @@name:	standalone.2f
! @@type:	F-free
! @@compilable:	yes
! @@linkable:	no
! @@expect:	success
SUBROUTINE STANDALONE_OK()
  INTEGER  A
  A = 1
  IF (A .NE. 0) THEN
    !$OMP FLUSH(A)
  ENDIF
  IF (A .NE. 0) THEN
    !$OMP BARRIER
  ENDIF
  IF (A .NE. 0) THEN
    !$OMP TASKWAIT
  ENDIF
  IF (A .NE. 0) THEN
    !$OMP TASKYIELD
  ENDIF
  GOTO 100
  100 CONTINUE
  !$OMP FLUSH(A)
  GOTO 200
  200 CONTINUE
  !$OMP BARRIER
  GOTO 300
  300 CONTINUE
  !$OMP TASKWAIT
  GOTO 400
  400 CONTINUE
  !$OMP TASKYIELD
END SUBROUTINE

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


program test_scalar_cholesky

  integer, parameter :: r64 = selected_real_kind (14),                         &
                        r32 = selected_real_kind (6)
  integer :: n

  integer :: i, j, k

  logical :: positive_definite

  real (r64), allocatable :: A (:,:), B (:,:), L (:,:)

  real start_time, stop_time

  read (*, *) n

  allocate ( A (1:n, 1:n), B (1:n, 1:n), L (1:n, 1:n) )

  ! -------------------------------------------------------------
  ! create a positive definite matrix A by setting A equal to the
  ! matrix-matrix product B B^T.  This normal equations matrix is
  ! positive-definite as long as B is full rank.
  ! -------------------------------------------------------------

    write (*,*) "Cholesky Factorization Tests"
    write (*,*) "   Matrix Dimension: ", n
    write (*,*) "   Block Size      : ", 1
    write (*,*) "   Indexing Base   : ", 1

   call random_number ( harvest = B )

  do i = 1, n
     do j = 1, n
        A (i,j) = 0.0
        do k = 1, n
           A (i,j) = A (i,j) + B (i,k) * B (j,k)
        end do
     end do
  end do

  ! factorization algorithms overwrite a copy of A, leaving
  ! the factor L in its place

  L = A

  call cpu_time ( start_time )
  call scalar_outer_product_cholesky ( n, L, positive_definite )
  call cpu_time ( stop_time )

  elapsed_time = stop_time - start_time
  write (*,*)  "  Cholesky Factorization time:   ", elapsed_time
  write (*,*)  "  Collective speed in megaflops: ",                            &
               ( (n**3) / 3.0 )  / ( 1.0E6 * elapsed_time )

  if ( positive_definite ) then
     call check_factorization ( n, A, L )
  else
     write (*,*) "factorization failed for non-positive semi-definite matrix"
  end if

end program test_scalar_cholesky

subroutine scalar_outer_product_cholesky ( n, A, positive_definite )

  integer, parameter :: r64 = selected_real_kind (14),                         &
                        r32 = selected_real_kind (6)
  integer    :: n
  real (r64) :: A (1:n, 1:n)
  logical    :: positive_definite

  integer    :: j, k

  positive_definite = .true.

  do j = 1, n

     if ( A (j,j) > 0.0 ) then

        A (j,j) = sqrt ( A (j,j) )

        A (j+1:n, j) = A (j+1:n, j) / A (j,j)

        do k = j+1, n
           A (k:n, k) = A (k:n, k) - A(k:n, j) * A (k, j)
        end do

     else

        positive_definite = .false.
        exit

     end if

  end do

end subroutine scalar_outer_product_cholesky

subroutine check_factorization ( n, A, L )

  integer, parameter :: r64 = selected_real_kind (14),                         &
                        r32 = selected_real_kind (6)
  integer    :: n

  real (r64) :: A (1:n, 1:n), L (1:n, 1:n)

  ! -----------------------------------------------------------------------
  ! Check the factorization by forming L L^T and comparing the result to A.
  ! The factorization is successful if the results are within the
  ! point-wise perturbation bounds of Demmel, given as Theorem 10.5 in
  ! Higham's "Accuracy and Stability of Numerical Algorithms, 2nd Ed."
  ! -----------------------------------------------------------------------

  real (r64) :: unit_roundoff, gamma_n1

  real (r64) :: d (1:n)

  real (r64) :: resid, max_ratio

  unit_roundoff = 2.0 ** (-53) ! IEEE 64 bit floating point,
  !  also epsilon (unit_roundoff) in Fortran

  gamma_n1      = (n * unit_roundoff) / (1.0 - n * unit_roundoff)

  max_ratio = 0.0

  do i = 1, n
     d (i) = sqrt ( A (i,i) )
  end do

  do i = 1, n
     do j = i+1, n

        resid = A (i,j)
        do k = 1, min (i,j)
           resid = resid - L (i,k) * L (j,k)
        end do
        resid = abs (resid)

        max_ratio = max ( max_ratio,                                           &
                          resid * (1 - gamma_n1) /                             &
                          ( gamma_n1 * d (i) * d (j) ) )

     end do
  end do

  if ( max_ratio <= 1.0 ) then
     write (*,*)   "factorization successful"
  else
     write (*,*) "factorization error exceeds bound by ratio:", max_ratio
  end if

end subroutine check_factorization


text:
The DeCarli Room is one of the meeting rooms on the second floor of the Memorial Union.  It is named in honor of Dean DeCarli.

## Image(DeCarli.jpg, thumbnail, 300, DeCarli as pictured on the Plaque)


text:
      MODULE coils
      USE coiltypes

! This module declares an allocatable array of coils for each type of coil.
! The actual allocation is done in the initialization routine.

         TYPE(modularcoil), DIMENSION(:), ALLOCATABLE :: modular
         TYPE(saddlecoil),  DIMENSION(:), ALLOCATABLE :: saddle
         TYPE(vfcoil),      DIMENSION(:), ALLOCATABLE :: vertical

      END MODULE coils

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
namespace Foo\Bar\User {

    use Foo\Bar\Command;
    use Foo\Bar\CommandId;
    use Foo\Bar\Name;
    use Foo\Bar\Age;

    command UserCommand : Command (CommandId) = RegisteredUser as registered-user { Name $name, ?Age $age}
        | UpdateUserName as update-user-name { Name $name }
        | DeleteUser;
}

ext: fpp

max_stars_repo_name: Mattin/fpp

max_issues_repo_name: Mattin/fpp

max_forks_repo_name: Mattin/fpp


text:
program cpu_time_01
use iso_fortran_env, only: dp=>real64
implicit none
real(dp) :: t1, t2
call cpu_time(t1)
print *, "Some computation"
call cpu_time(t2)
print *, "Total time: ", t2-t1
end program

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


!=======================================================================
!     RECIPROCAL
!=======================================================================

module reciprocal_inp
  !     k-space variables :
  use controls  !KJ 8/06
  use struct, nphstr => nph
  use kklist,only: nkp,usesym,nkx,nky,nkz,ktype
  use strfacs,only: streta,strrmax,strgmax,init_strfacs
  implicit none
  integer icorehole
  real*8 streimag ! additional broadening for calculation KKR structure factors ; not recommended
  character(*),parameter,private :: filename='reciprocal.inp'

contains

  subroutine reciprocal_write
    !KJ next file added 8/06
    integer i
    open (file=filename, unit=3, status='unknown')
    !       in which space are we?
    write(3,10) 'spacy'
    write(3,20) spacy
    if(spacy.eq.0) then
       write(3,10) 'lattice vectors  (in A, in Carthesian coordinates)'
       write(3,30) a1
       write(3,30) a2
       write(3,30) a3
       write(3,10) 'Volume scaling factor (A^3); eimag; core hole'
       write(3,30) dble(-1),dble(0),dble(1)
       write(3,10) 'lattice type  (P,I,F,R,B,CXY,CYZ,CXZ)'
       write(3,10) latticename
       write(3,10) '#atoms in unit cell ; position absorber ; corehole?'
       write(3,20) nats,absorber,icorehole
       write(3,10) '# k-points total/x/y/z ; ktype; use symmetry?'
       write(3,*) nkp,nkx,nky,nkz,ktype,usesym  ! format line 20 limits integer to 4 positions - not enough for nkp!
       write(3,10) 'ppos'
       do i=1,nats
          write(3,30) ppos(:,i)
       enddo
       write(3,10) 'ppot'
       !KJ bugfix 5/2012: It's important not to use formatting when there are more atoms than fit on one line!!
       write(3,*) ppot
       write(3,10) 'streta,strgmax,strrmax'
       write(3,30) streta,strgmax,strrmax
    endif
    close(3)
    ! standard formats for string, integers and real numbers
10  format(a)
20  format (20i4)
30  format (9f13.5)
  end subroutine reciprocal_write

  subroutine reciprocal_read(celvin)
    use struct, nphstr => nph
    integer i
    real*8,intent(out) :: celvin
    open (3,file=filename,status='unknown',err=167)
    read(3,*,end=167,err=167)
    read(3,*,end=167,err=167) spacy
    if(spacy.eq.0) then
       read(3,*) ; read(3,*) a1(:)
       read(3,*) a2(:)
       read(3,*) a3(:)
       read(3,*) ; read(3,*) celvin,streimag,cholestrength
       read(3,*) ; read(3,*) latticename
       lattice=latticename(1:1)
       read(3,*) ; read(3,*) nats,absorber,icorehole
       read(3,*) ; read(3,*) nkp,nkx,nky,nkz,ktype,usesym
       read(3,*)
       !Careful: the next statement used to be "if size(ppot).eq.0".  However, on ifort size(ppot)=0 but on gfortran it =1!!
       !Hence the new instruction.
       !I wish if(allocated(ppot)) would work here; I don't understand why it doesn't.
       if(size(ppot).lt.nats) call init_struct(nats) !KJ 7-09 bugfix call this only once ; I can't seem to use "allocated(ppos)" here?
       do i=1,nats
          read(3,*) ppos(:,i)
       enddo
       read(3,*) ; read(3,*) ppot
       read(3,*) ; read(3,*) streta,strgmax,strrmax
       if(icorehole.eq.1) then
          corehole=.true.
       else
          corehole=.false.
       endif
    endif
    return
167 spacy=1
    return
  end subroutine reciprocal_read

  subroutine reciprocal_init
    call init_controls
    call init_strfacs
    icorehole = 1  ! use core hole
    streimag = dble(0) ! no extra broadening for KKR struc factors
    cholestrength = dble(1) ! don't mess with core hole
  end subroutine reciprocal_init

end module reciprocal_inp