      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),CNV,STRX,STRTOT
      INTEGER*2 IPTR(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,IFLG,
     :        I0,I1,NT,ISTATE,J,K,N,NINDX,NNINDX,KPOSU,KPOSL,IERR
      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
      DO 2 I=1,NDX
   2      IPTR(I)=I+1
      IPTR(NDX)=0
      KOFF(1)=0
      DO 3 I=1,NDT
   3      KOFF(I+1)=KOFF(I)+I
      ISTATE=-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)',IOSTAT=IERR)
     :               FRQX,STRX,IQNFMT,IQNUX,IQNLX,INDX
      IF(IERR.NE.0) STOP 'bad STR file' 
      IF(INDX.LE.0) INDX=1
      NINDX=0
      NNINDX=0
      NT=0
      STRTOT=0.
      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=1+IFHALF
C  begin master loop
  10  IQNU=IQNUX
      IQNL=IQNLX
      FRQC=FRQX
      STR=STRX
      IF(ISTATE.EQ.0) THEN
          READ(LUSTR,'(2E15.5,6X,A,A,I5)',IOSTAT=IERR)
     :            FRQX,STRX,IQNUX,IQNLX,INDX
          IF(INDX.LE.0) INDX=1
          IF(IERR.NE.0) THEN
              ISTATE=1
          ENDIF
      ENDIF
      IF(ISTATE.LE.0) THEN
        IF(IQNUX.EQ.IQNU.AND.IQNLX.EQ.IQNL) THEN
C        same quanta, different dipole
          IF(INDX.GT.NINDX) THEN
              IF(INDX.GT.NDT) THEN
                  WRITE(*,*) 'NINDX, NDT =',INDX,NDT
                  STOP
              ENDIF
              NINDX=INDX
              I=NNINDX
              NNINDX=KOFF(INDX+1)
              N=NNINDX-I
              I=I+1
C              zero out extra A and B values
              DO 11 J=1,2
                 KPTR=HEAD(J)
C                while ..
  12             IF(KPTR.GT.0) THEN
                    K=(KPTR-1)*NDTT+I
                    CALL DCOPY(N,ZERO,I0,A(K),I1)
                    CALL DCOPY(N,ZERO,I0,B(K),I1)
                    KPTR=IPTR(KPTR)
                    GO TO 12
                 ENDIF
  11          CONTINUE
          ENDIF
          NT=NT+1
          STRV(NT)=STRX
          INDXV(NT)=INDX
          STRTOT=STRTOT+STRX
          ISTATE=0
          GO TO 10
        ENDIF
      ENDIF
      IF(ISTATE.NE.2) THEN
          IQFU=IVAL(IQNU(NQN:NQN+1))
          IQFL=IVAL(IQNL(NQN:NQN+1))
          IQFX=MAX(IQFL,IQFU)
          IF(IQFX.LT.IQFMX) THEN
              ISTATE=2
              IQFX=IQFMX+1
          ENDIF
      ENDIF
C     while .. purge data for old F
  15  IF(IQFX.GT.IQFMX) THEN
          KPTR=HEAD(1)
C         while .. more data for F
  20      IF(KPTR.GT.0) THEN
              K=(KPTR-1)*NDTT
              DO 25 I=1,NINDX
                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
          IF(HEAD(1).EQ.0) IQFMX=IQFX
          IF(IQFX.EQ.IQFMX) WRITE(*,*) 'Starting F=',IQFX
          GO TO 15
      ENDIF
      IF(ISTATE.EQ.2) STOP
      STRTOT=ABS(STRTOT)*CNV
      IF(FRQC.LT.BIG*STRTOT) 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
          IFLG=IQFU-IQFMX
          IPOSU=LPOSN(IQNU,IFLG,HEAD,IPTR,IQN)
          KPOSU=(IPOSU-1)*NDTT
          IF(IFLG.GT.0) THEN
              CALL DCOPY(NNINDX,ZERO,I0,A(KPOSU+1),I1)
              CALL DCOPY(NNINDX,ZERO,I0,B(KPOSU+1),I1)
          ENDIF
          IFLG=IQFL-IQFMX  
          IPOSL=LPOSN(IQNL,IFLG,HEAD,IPTR,IQN)
          KPOSL=(IPOSL-1)*NDTT
          IF(IFLG.GT.0) THEN
              CALL DCOPY(NNINDX,ZERO,I0,A(KPOSL+I),I1)
              CALL DCOPY(NNINDX,ZERO,I0,B(KPOSL+I),I1)
          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)
              TA=TA/FRQC
              TB=TB/FRQC
              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
      STRV(1)=STRX
      INDXV(1)=INDX
      STRTOT=STRX
      IF(ISTATE.EQ.1) THEN
          ISTATE=2
          IQFX=IQFMX+1
      ENDIF
      GO TO 10
      END
      INTEGER FUNCTION LPOSN(IQNU,IDX,HEAD,IPTR,IQN)
      INTEGER IDX,HEAD(3),KPTR
      INTEGER*2 IPTR(*)
      CHARACTER*12 IQN(*),IQNU
      IF(IDX.NE.0) IDX=1
      IDX=2-IDX
      KPTR=HEAD(IDX)
 10   IF(KPTR.GT.0) THEN
          LPOSN=KPTR
          KPTR=IPTR(KPTR)
          IF(IQN(LPOSN).NE.IQNU) GO TO 10
          IDX=0
          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(IDX)
      HEAD(IDX)=LPOSN
      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

