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************************************************************************
      INTEGER NDX,NDT,NDTT
      PARAMETER (NDX=4000,NDT=16,NDTT=(NDT*NDT+NDT)/2)
      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
      DATA CEXT/'str ','stk '/
      BIG=100.
	EPS=1.E-12
      CNV=0.503403
      ZERO(1)=0.
      I0=0
      I1=1
      LUSTR=11
      LUSTK=13
      HEAD(1)=0
      HEAD(2)=0
      HEAD(3)=1
      NTP=1
      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
      CALL FILGET(CFIL,CEXT,2)
      CALL OPENLU(LUSTR,CFIL(1),'R')
      CALL OPENLU(LUSTK,CFIL(2),'W')
      READ(LUSTR,'(2E15.5,I5,1X,A,A,I5)',END=50)
     :               FRQX,STRX,IQNFMT,IQNUX,IQNLX,INDX
      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  begin master loop
  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(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     while ..
  15  IF(IQFX.GT.IQFMX) THEN
          WRITE(*,*) 'Starting F=',IQFX
          KPTR=HEAD(1)
C         while ..
  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
                     WRITE(LUSTK,'(2E15.6,1X,A,2I5)')
     &                     TA,TB,IQN(KPTR),J,I
                  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) STOP
      STR=ABS(STR)*CNV
      IF(FRQC.LT.BIG*STR) THEN
C        degenerate
          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))
              WRITE(LUSTK,'(F15.4,E15.6,1X,A,A,I5)') 
     &                       FRQC,STR,IQNU,IQNL,I
            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
      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
  50  DONE=0
      GO TO 12
      END
      INTEGER FUNCTION LPOSN(IQNU,IDX,HEAD,IPTR,NAB,IQN)
      INTEGER IDX,HEAD(3),KPTR,IDXX
      INTEGER*2 IPTR(*),NAB(*)
      CHARACTER*12 IQN(*),IQNU
      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
      RETURN
      END
      SUBROUTINE STARKS(A,B,IQFX,IQDEL,IFHALF)
      INTEGER IQFX,IQDEL,IFHALF,IDGN
      REAL*8 A,B,F,FSQ
      IDGN=IQFX+IQFX-IFHALF
      F=0.5*IDGN
      FSQ=F*F
      IF(IQDEL.EQ.0) THEN
          B=A/((F+FSQ)*(IDGN+1))
          A=0
      ELSE
          A=A/(4.*F-1./F)
          B=-A/FSQ
      ENDIF
      RETURN
      END
      INTEGER FUNCTION IVAL(STR)
      CHARACTER*2 STR
	INTEGER I
	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')
      RETURN
      END