c$DEBUG
C----------------------------------------------------------------------
C  Downloaded from http://spec.jpl.nasa.gov/ftp/pub/calpgm on 10.XI.98
c
c  MOdifications to original (for executable lines marked by zk in
c  the comment column):
c
c  - Explicit inclusion of necessary routines from SLIB.F and BLAS.F
c  - Debugging of second order calc.
c  - Touching up for better readability of output and of listing (with
c    apologies for the 'deglev' hack)
c
C    v.27.XI.1998                                 Zbigniew KISIEL
c                                                 kisiel@ifpan.edu.pl
C----------------------------------------------------------------------
C
      PROGRAM STARK
C
C   Copyright (C) 1989, California Institute of Technology
C   All rights reserved.  U. S. Government Sponsorship under
C   NASA Contract NAS7-918 is acknowledged.
C
C   Herbert M. Pickett, 19 December 1989
C
C     get command line file names (.str = calc input, .stk = output)
C     open .str and read first line for NQN
C     DO UNTIL end of file on .str
C        IF Fup or Flow > Fmax THEN
C            printout stark coefficients for Fmax-2
C        ENDIF
C        find lower level or create new entry
C        find upper level or create new entry
C        IF degeneracy THEN
C            print info
C        ELSE
C            calc stark effect and add to upper and lower
C        ENDIF
C        read in line from .str
C     END DO
C************************************************************************
c
      INTEGER NDX,NDT,NDTT
c...HMP dimensioning:
c     PARAMETER (NDX=4000,NDT=16,NDTT=(NDT*NDT+NDT)/2)
c...custom dimensioning
      PARAMETER (NDX=4000,NDT=64,NDTT=(NDT*NDT+NDT)/2)
c
      CHARACTER CFIL(2)*64, CEXT(2)*4, IQN(NDX)*12, IQNL*12, IQNU*12,
     :           IQNLX*12, IQNUX*12
      REAL*8 A(NDTT*NDX),B(NDTT*NDX),TA,TB,FRQC,FRQX,STR,EPS,
     :       STRV(NDT),BIG,ZERO(2),STMP,CNV,STRX
      INTEGER*2 IPTR(NDX),NAB(NDX),INDXV(NDT),KOFF(NDT+1)
      INTEGER IQNFMT,IFHALF,IQFU,IQFL,HEAD(3),KPTR,IQFMX,IQFX
     :        ,LPOSN,IPOSU,IPOSL,I,LUSTR,LUSTK,NQN,IVAL,INDX,NTP
     :        ,I0,I1,NT,DONE,J,K,N,INDXMX,KPOSU,KPOSL
      character deglev(ndtt)*12,dflag*1                                  zk
c
      DATA CEXT/'str ','stk '/
c
c...BIG controls picking out of near degeneracy: set this as necessary
      BIG=10000.
c
      EPS=1.E-12
c
c...CNV is the conversion factor to use field in V/cm, dipoles in Debyes and
c   frequency in MHz (=10/(2.99792458*6.6260755); 1986 constants, only the
c   number from the Planck constant is subject to change)
c   Note that the Debye was defined in the CGS era as 10^(-10) of the
c   dipole moment from two charges equal to +- 1 cgs charge unit (franklin)
c   separated by 1 Angstrom = 10^(-18) franklin cm;
c   1 Franklin = (1/c(CGS)) Coulomb = 3.33564095*10^(-10) C so that
c           1D = 3.33564095*10^(-30) C
c
      CNV=0.50341125D0
c
      ZERO(1)=0.
      I0=0
      I1=1
      LUSTR=11
      LUSTK=13
      HEAD(1)=0
      HEAD(2)=0
      HEAD(3)=1
      NTP=1
      ndglev=0                                                           zk
c
      DO 2 I=1,NDX
   2      IPTR(I)=I+1
      KOFF(1)=0
      DO 5 I=1,NDT
   5      KOFF(I+1)=KOFF(I)+I
      IPTR(NDX)=0
      IQFMX=1
      NQN=1
      DONE=-1
c
      CALL FILGET(CFIL,CEXT,2)
      CALL OPENLU(LUSTR,CFIL(1),'R')
      CALL OPENLU(LUSTK,CFIL(2),'W')
c
c...Block of additional descriptive output (zk)
c
      write(lustk,5550)
5550  format(' Pure second order Stark effect:'/
     * '      E.level = E.o + ( A.g + B.g M^2 ) mu.g^2 E^2'/
     * ' Mixed order:'/
     * '      E.level = (E1+E2)/2  +- [ (DE/2)^2 ',
     * '+ mu.g^2 E^2 |xi|^2 M^2 ]^(1/2)'//
     * '              = E1-DE/2+[...]^(1/2) and'/
     * '                E2+DE/2-[...]^(1/2) for E1>E2, DE=E1-E2'//
     * ' * E1,E2,DE should in general be inclusive of second order '
     * 'Stark contributions'/
     * ' * Asterisk denotes levels implicated in near degeneracy',
     * ' for which'/'   A and B are evaluated without the ',
     * 'connecting (|Phi.Zg|) term'/
     * ' * Sum field contributions for more ',
     * 'than one dipole component'/
     * ' * Dipoles in D, energies in MHz, field in V/cm'/
     * ' * In the .INT file set mu.g to 1 (only one at a time)',
     * ' and STRFLG=1'/)
      do 5554 i=64,1,-1
        if(cfil(1)(i:i).eq.'.')goto 5555
5554  continue
5555  cfil(1)=cfil(1)(1:i-1)//'.int'
      open(3,file=cfil(1),status='old',err=5552)
      read(3,'(a)',err=5552)cfil(1)
      write(lustk,5553)cfil(1)
5553  format(1x,72(1h-)/1x,a/1x,72(1h-))
      close(3)
5552  write(lustk,5551)
5551  format(1x//'           A.g           B.g       level',25x,'g'/)
c
c
c...read the first transition and work out the position of the F quantum number
c   in fields IQNUX and IQNLX by using the formatting information in IQNFMT
c   The first character of F is defined by NQN
c
c   Note that STRX is the reduced transition dipole matrix element equal to
c   square root of the linestrength
c
      READ(LUSTR,'(2E15.5,I5,1X,A,A,I5)',END=50)
     :               FRQX,STRX,IQNFMT,IQNUX,IQNLX,INDX
      if(frqx.eq.0.0)frqx=1.d-05                                         zk
      strx=strx/dsqrt(frqx)                                              zk
      IF(INDX.LE.0) INDX=1
      NT=1
      STRV(1)=STRX
        INDXV(1)=INDX
        INDXMX=INDX
      IFHALF=IQNFMT/10
      NQN=IQNFMT-IFHALF*10
      IF( MOD(IFHALF/10,5) .EQ.NQN) NQN=1
      NQN=NQN+NQN-1
      IFHALF=MOD(IFHALF,2)
      IQFMX=IQFMX+IFHALF
c
C  begin master loop
c
  10  IQNU=IQNUX
      IQNL=IQNLX
      FRQC=FRQX
      STR=STRX
      READ(LUSTR,'(2E15.5,6X,A,A,I5)',END=50)
     :            FRQX,STMP,IQNUX,IQNLX,INDX
      if(frqx.eq.0.0)frqx=1.d-05                                         zk
      stmp=stmp/dsqrt(frqx)                                              zk
      IF(INDX.LE.0) INDX=1
      IF(IQNUX.EQ.IQNU.AND.IQNLX.EQ.IQNL) THEN
          NT=NT+1
          INDXV(NT)=INDX
          IF(INDXMX.LT.INDX) INDXMX=INDX
          STRV(NT)=STMP
          STRX=STRX+STMP
          GO TO 10
      ENDIF
  12  IQFU=IVAL(IQNU(NQN:NQN+1))
      IQFL=IVAL(IQNL(NQN:NQN+1))
      IQFX=MAX(IQFL,IQFU)
c
C     while ..
c
  15  IF(IQFX.GT.IQFMX) THEN
          WRITE(*,*) 'Starting F=',IQFX
          write(lustk,'(1x)')                                            zk
          KPTR=HEAD(1)
c
C         while ..
c
  20      IF(KPTR.GT.0) THEN
              N=NAB(KPTR)
              K=(KPTR-1)*NDTT
              DO 25 I=1,N
                DO 27 J=1,I
                    K=K+1
                  TA=A(K)
                  TB=B(K)
                  IF(ABS(TA)+ABS(TB).GT.EPS) THEN
                     dflag=' '                                           zk
                     if(ndglev.ne.0)then                                 zk
                       mzk=ndglev                                        zk
                       do 5560 izk=1,mzk                                 zk
                         if(iqn(kptr).eq.deglev(izk))then                zk
                           dflag='*'                                     zk
                           ndglev=ndglev-1                               zk
                           do 5562 nzk=izk,ndglev                        zk
                             deglev(nzk)=deglev(nzk+1)                   zk
5562                       continue                                      zk
                           goto 5561                                     zk
                         endif                                           zk
5560                   continue                                          zk
5561                   continue                                          zk
                     endif                                               zk
                     WRITE(LUSTK,
     *                   '(1x,a1,1P(E14.6,1x,e16.6),1X,A,15x,2I5)')      zk
     &                     dflag,TA,TB,IQN(KPTR),J,I                     zk
                  ENDIF
  27              CONTINUE
  25            CONTINUE
              I=IPTR(KPTR)
              IPTR(KPTR)=HEAD(3)
              HEAD(3)=KPTR
              KPTR=I
              GO TO 20
          ENDIF
          IQFMX=IQFMX+1
          HEAD(1)=HEAD(2)
          HEAD(2)=0
          GO TO 15
      ENDIF
      IF(DONE.GT.0) then
        write(lustk,'(1x/2(1x,72(1h-)/))')                               zk
        STOP
      endif                                                              zk
      STR=ABS(STR)*CNV
      IF(FRQC.LT.BIG*STR) THEN
c
C        degenerate
c
          DO 30 I=1,NT
            TA=STRV(I)*CNV
            TA=TA*TA
            IF(TA.GT.EPS) THEN
              CALL STARKS(TA,TB,IQFX,IQFU-IQFL,IFHALF)
              STR=SQRT(MAX(TA,TB))
              str=str*dsqrt(frqc)                                        zk
              WRITE(LUSTK,'('' DE='',f10.4,'' |xi|='',1PE13.6            zk
     *                    ,1X,A,''<->'',A,I5)')                          zk
     &                       FRQC,STR,IQNU,IQNL,I
              deglev(ndglev+1)=iqnu                                      zk
              deglev(ndglev+2)=iqnl                                      zk
              ndglev=ndglev+2                                            zk
            ENDIF
 30         CONTINUE
      ELSE
          N=KOFF(INDXMX+1)
          IPOSU=LPOSN(IQNU,IQFMX-IQFU,HEAD,IPTR,NAB,IQN)
          KPOSU=(IPOSU-1)*NDTT
          IF(NAB(IPOSU).LT.INDXMX) THEN
            I=NAB(IPOSU)
            IF(I.GT.0) I=KOFF(I+1)
            J=N-I
            I=I+1
            CALL DCOPY(J,ZERO,I0,A(KPOSU+I),I1)
            CALL DCOPY(J,ZERO,I0,B(KPOSU+I),I1)
            NAB(IPOSU)=INDXMX
          ENDIF
          IPOSL=LPOSN(IQNL,IQFMX-IQFL,HEAD,IPTR,NAB,IQN)
          KPOSL=(IPOSL-1)*NDTT
          IF(NAB(IPOSL).LT.INDXMX) THEN
            I=NAB(IPOSL)
            IF(I.GT.0) I=KOFF(I+1)
            J=N-I
            I=I+1
            CALL DCOPY(J,ZERO,I0,A(KPOSL+I),I1)
            CALL DCOPY(J,ZERO,I0,B(KPOSL+I),I1)
            NAB(IPOSL)=INDXMX
          ENDIF
          DO 35 I=1,NT
            DO 40 J=1,I
              TA=STRV(I)*CNV
              IF(I.EQ.J) THEN
                TA=TA*TA
              ELSE
                TA=2.*TA*STRV(J)*CNV
              ENDIF
              CALL STARKS(TA,TB,IQFX,IQFU-IQFL,IFHALF)
              K=INDXV(J)+KOFF(INDXV(I))
              A(K+KPOSU)=A(K+KPOSU)+TA
              A(K+KPOSL)=A(K+KPOSL)-TA
              B(K+KPOSU)=B(K+KPOSU)+TB
              B(K+KPOSL)=B(K+KPOSL)-TB
  40        CONTINUE
  35      CONTINUE
      ENDIF
c
      NT=1
      INDXV(1)=INDX
      INDXMX=INDX
      STRV(INDX)=STMP
      STRX=STMP
      IF(DONE.LT.0) GO TO 10
      DONE=1
      IQFX=IQFMX+1
      GO TO 15
C
c...end of .STR file, return to 12 to complete calculations
c
  50  DONE=0
      GO TO 12
c
      END
c
c--------------------------------------------------------------------------
c
      INTEGER FUNCTION LPOSN(IQNU,IDX,HEAD,IPTR,NAB,IQN)
c
      INTEGER IDX,HEAD(3),KPTR,IDXX
      INTEGER*2 IPTR(*),NAB(*)
      CHARACTER*12 IQN(*),IQNU
c
      IDXX=2-IDX
      KPTR=HEAD(IDXX)
 10   IF(KPTR.GT.0) THEN
          LPOSN=KPTR
          KPTR=IPTR(KPTR)
          IF(IQN(LPOSN).NE.IQNU) GO TO 10
          RETURN
      ENDIF
C  take from top of free and put at top of new
      LPOSN=HEAD(3)
      IF(LPOSN.EQ.0) STOP 'LPOSN=0'
      HEAD(3)=IPTR(LPOSN)
      IPTR(LPOSN)=HEAD(IDXX)
      HEAD(IDXX)=LPOSN
      NAB(LPOSN)=0
      IQN(LPOSN)=IQNU
c
      RETURN
      END
c
c--------------------------------------------------------------------------
c
      SUBROUTINE STARKS(A,B,IQFX,IQDEL,IFHALF)
c
c...IQFX is the current value of F which is greater or equal to that
c   for the other level: note that with such F calculation of P and R dipole
c   elements using G&C eq.10.68,71-73 becomes symmetric
c   IQDEL is DeltaF between the two levels
c   IFHALF =0,1 for even,odd number of 1/2 in F
c
      INTEGER IQFX,IQDEL,IFHALF,IDGN
      REAL*8 A,B,F,FSQ
c
      IDGN=IQFX+IQFX-IFHALF
      F=0.5d0*real(IDGN)
      FSQ=F*F
      IF(IQDEL.EQ.0) THEN
          B=A/((F+FSQ)*real(IDGN+1))
          A=0.0d0
      else
          A=A/(4.d0*F-1.d0/F)
c                              equivalent to A=A*F/(4.d0*F^2-1.d0)       zk
          B=-A/FSQ
      ENDIF
c
      RETURN
      END
c
c--------------------------------------------------------------------------
c
      INTEGER FUNCTION IVAL(STR)
c
c...calculate the decimal value corresponding to the two numeric characters
c   in STR, with use of 'A' and above for 11,12 etc. in first digit
c
      CHARACTER*2 STR
      INTEGER I
c
      I=ICHAR(STR(1:1))
      IVAL=I-ICHAR('0')
      IF(IVAL.GT.9) IVAL=10+I-ICHAR('A')
      IF(IVAL.LT.0) IVAL=0
      IF(STR(2:2).NE.' ') IVAL=10*IVAL+ICHAR(STR(2:2))-ICHAR('0')
c
      RETURN
      END
c
c--------------------------------------------------------------------------
c
      SUBROUTINE FILGET(CFIL,CEXT,NFILE)
C
C     Get NFILE names with extent given in CEXT
C     Put names in CFIL
C
C
      CHARACTER*(*) CFIL(*),CEXT(*)
      INTEGER NFILE
      CHARACTER EXT*3,EXTQ*3,FIL*63,EXTDOT*1
      INTEGER NLEN,J,K,IEXT,ICHAR,IEXTQ
C     define character which delimits EXT field
      EXTDOT='.'
C  get first parameter
      WRITE(*,*) ' Enter file name '
      READ(*,100) FIL
      NLEN=INDEX(FIL,EXTDOT)
      IF(NLEN.EQ.0) THEN
C         extdot not found
          NLEN=INDEX(FIL,' ')
          IF(NLEN.LE.1) STOP
          FIL(NLEN:NLEN)=EXTDOT
      ENDIF
C     set up defaults
      DO 10 J=1,NFILE
          EXT=CEXT(J)
 10       CFIL(J)=FIL(1:NLEN)//EXT
C     find non default values
      DO 20 I=2,20
          WRITE(*,*) ' Enter file name '
          READ(*,100,END=50) FIL
          K=INDEX(FIL,EXTDOT)
          IF(K.EQ.0) GO TO  50
          EXT=FIL(K+1:)
          DO 30 J=1,NFILE
              EXTQ=CEXT(J)
              DO 35 K=1,3
                  IEXT=ICHAR(EXT(K:K))
                  IEXTQ=ICHAR(EXTQ(K:K))
                  IF(IEXT.EQ.IEXTQ) GO TO 35
                  IF(IEXT+32.EQ.IEXTQ) GO TO 35
                      GO TO 30
35                CONTINUE
              CFIL(J)=FIL
 30           CONTINUE
 20       CONTINUE
 50   WRITE(*,'(1X,A)') (CFIL(J),J=1,NFILE)
      RETURN
 100  FORMAT(A)
c
      END
c
c--------------------------------------------------------------------------
c
      SUBROUTINE OPENLU(LU,FNAME,CHSW)
c
C     fancy file opening and closing
c
      INTEGER LU,LUORG,LUIN,IREC,IOCHK
      PARAMETER (LUORG=10)
      CHARACTER*8 LULIST
      CHARACTER*(*) FNAME,CHSW
      CHARACTER*64  FIL
      SAVE LULIST
      DATA LULIST /'CCCCCCCC'/
      LUIN=LU
      IF(CHSW(1:1).EQ.'C') THEN
          LU=LU-LUORG
          IF(LU.GT.0.AND.LU.LE.LEN(LULIST)) THEN
              IF(LULIST(LU:LU).EQ.'R') THEN
                  REWIND (LUIN)
              ELSE IF(LULIST(LU:LU).EQ.'W') THEN
                  ENDFILE (LUIN)
              ENDIF
              CLOSE(LUIN,ERR=50)
              LULIST(LU:LU)='C'
          ENDIF
      ELSE
          LU=INDEX(LULIST,'C')
          IF(LU.EQ.0) THEN
              WRITE(*,*) ' Too many open files'
              STOP
          ENDIF
          LULIST(LU:LU)=CHSW(1:1)
          LU=LU+LUORG
          IF(CHSW.EQ.'R') THEN
            OPEN(LU,FILE=FNAME,ERR=99,STATUS='OLD',IOSTAT=IOCHK)
          ELSE IF(CHSW.EQ.'RL') THEN
            OPEN(LU,FILE=FNAME,ERR=80,STATUS='OLD')
          ELSE IF(CHSW.EQ.'RB') THEN
            OPEN(LU,FILE=FNAME,ERR=85,STATUS='OLD',
     :           FORM='UNFORMATTED')
          ELSE IF(CHSW.EQ.'W') THEN
            OPEN(LU,FILE=FNAME,ERR=99,STATUS='UNKNOWN',IOSTAT=IOCHK)
          ELSE IF(CHSW.EQ.'WP') THEN
c
C           allow others to read my partially written file
c
            OPEN(LU,FILE=FNAME,ERR=99,STATUS='UNKNOWN',IOSTAT=IOCHK)
          ELSE IF(CHSW.EQ.'WB') THEN
            OPEN(LU,FILE=FNAME,ERR=85,STATUS='UNKNOWN',
     :           FORM='UNFORMATTED')
          ELSE IF(CHSW.EQ.'S') THEN
            IREC=LUIN*8
            OPEN(LU,STATUS='SCRATCH',ACCESS='DIRECT',RECL=IREC,ERR=99
     :        ,IOSTAT=IOCHK)
          ELSE
            WRITE(*,*) ' OPENLU unknown option =',CHSW
            STOP
          ENDIF
      ENDIF
  50  RETURN
c
C    try to open from default directory for parameter labels
c
  80  FIL='/SPECTRA/'//FNAME
      OPEN(LU,FILE=FIL,ERR=90,STATUS='OLD')
      RETURN
  85  LU=0
      RETURN
  90  WRITE(*,*) ' Trouble opening ',FIL
      RETURN
  99  WRITE(*,*) ' Trouble opening ',FNAME
      WRITE(*,*) ' IOSTAT = ',IOCHK
c
      STOP
      END
c
c--------------------------------------------------------------------------
c
      SUBROUTINE  DCOPY(N,SX,INCX,SY,INCY)
C
C     COPIES A VECTOR, X, TO A VECTOR, Y.
C
      REAL*8  SX(0:*),SY(0:*),SA
      INTEGER I,INCX,INCY,N,M
      INTEGER*4 IX,IY
C
      M=N-1
      IF(M.LT.0)RETURN
      IF(INCX.EQ.0) THEN
        SA=SX(0)
        IF(INCY.NE.1) THEN
C
C        CODE FOR FILL AND INCREMENT NOT EQUAL TO 1
C
          IY = 0
          IF(INCY.LT.0) IY = -M*INCY
          DO 10 I = 0,M
            SY(IY) = SA
            IY = IY + INCY
   10       CONTINUE
        ELSE
C
C        CODE FOR FILL AND INCREMENT EQUAL TO 1
C
C
          DO 20 I = 0,M
            SY(I) = SA
   20       CONTINUE
        ENDIF
      ELSE IF(INCX.EQ.INCY) THEN
        IF(INCY.EQ.1) THEN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
          DO 30 I = 0,M
            SY(I) = SX(I)
   30       CONTINUE
C
C        CODE FOR BOTH INCREMENTS EQUAL
C
        ELSE
          IX = 0
          IF(INCX.LT.0)IX = -M*INCX
          DO 40 I = 0,M
            SY(IX) = SX(IX)
            IX = IX + INCX
   40     CONTINUE
        ENDIF
      ELSE
C
C        CODE FOR UNEQUAL INCREMENTS
C
        IX = 0
        IF(INCX.LT.0)IX = -M*INCX
        IY = 0
        IF(INCY.LT.0)IY = -M*INCY
        DO 50 I = 0,M
          SY(IY) = SX(IX)
          IX = IX + INCX
          IY = IY + INCY
   50     CONTINUE
      ENDIF
      RETURN
      END
      FUNCTION DDOT(N,SX,INCX,SY,INCY)
C
C     FORMS THE DOT PRODUCT OF TWO VECTORS.
C
      REAL*8  DDOT,SX(0:*),SY(0:*)
      INTEGER I,INCX,INCY,N,M
      INTEGER*4 IX,IY
C
      DDOT = 0
      M=N-1
      IF(M.LT.0)RETURN
      IF(INCX.EQ.INCY) THEN
        IF(INCY.EQ.1) THEN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
          DO 10 I = 0,M
            DDOT = DDOT + SX(I)*SY(I)
   10       CONTINUE
C
C        CODE FOR BOTH INCREMENTS EQUAL
C
        ELSE
          IX = 0
          IF(INCX.LT.0) IX = -M*INCX
          DO 20 I = 0,M
            DDOT = DDOT + SX(IX)*SY(IX)
            IX = IX + INCX
   20       CONTINUE
        ENDIF
      ELSE
C
C        CODE FOR UNEQUAL INCREMENTS
C
        IX = 0
        IF(INCX.LT.0) IX = -M*INCX
        IY = 0
        IF(INCY.LT.0) IY = -M*INCY
        DO 30 I = 0,M
          DDOT = DDOT + SX(IX)*SY(IY)
          IX = IX + INCX
          IY = IY + INCY
   30     CONTINUE
      ENDIF
c
      RETURN
      END
c
c--------------------------------------------------------------------------
c--------------------------------------------------------------------------