C N. R. BADNELL  XTRCTG            UoS v1.3                     29/10/14
C
      PROGRAM XTRCTG
      IMPLICIT REAL*8 (A-H,O-Z)
      LOGICAL BFORM
C
C BASIC USAGE:
C
C READ AN (R-MATRIX COLLISION STRENGTH) XOUT FILE PRODUCED BY XTRCT.F
C AND CONVOLUTE IT WITH AN EWIDTH FWHM GAUSSIAN BETWEEN EMIN AND EMAX
C USING NGAUSS POINTS.
C
C THE USER INPUT ENERGIES ARE IN UNITS (DEFAULT =1, RYD)
C THE OUTPUT TO XOUTG IS THEN IN UNITS ALSO.
C (THE UNITS IN XOUT ARE ASSUMED TO BE RYDBERGS STILL!)
C
C IF THE USER INPUTS AN INITIAL STATISTICAL WEIGHT IWT, THEN THE CODE
C COVERTS TO, CONVOLUTES AND OUTPUTS CROSS SECTIONS INSTEAD.
C
C ADVANCED OPTIONS:
C
C ILOG.LT.0 USES A LOGARITHMIC ENERGY MESH STARTING AT 10**ILOG RYD.
C TPAR.GT.0 USES STORAGE RING (COOLER) GAUSSIAN OF (FWHM) WIDTH
C           4*SQRT(LOG(2)*TPAR)*SQRT(E) WHERE TPAR IS THE PARALLEL TEMP.
C
      PARAMETER (ONE=1)
      PARAMETER (TWO=2)
      PARAMETER (TEN=10)
      PARAMETER (EINF=1000000)
C
      PARAMETER (NMXE=100000)
C
      DIMENSION EVEC(NMXE),VEC(NMXE)
C
      NAMELIST/GAUSS/NGAUSS,ILOG,E0,EMIN,EMAX,EWIDTH,TPAR,IWT,UNITS
C
      PIE=ACOS(-ONE)
      CONG=2*SQRT(LOG(TWO))
      CONX=28.00284D0                                        !A0**2 (MB)
      CONX=CONX*PIE                                       !PI*A0**2 (MB)
C
      INQUIRE (FILE='dxtrctg',EXIST=BFORM)
      IF(BFORM)THEN
        OPEN(5,FILE='dxtrctg',STATUS='UNKNOWN')      !CAN REDIRECT UNIT5
      ENDIF
C
      UNITS=1                                                  !RYDBERGS
      E0=0
      EMIN=0
      EMAX=-1
      EWIDTH=-1
      TPAR=-1
      NGAUSS=0
      ILOG=1                                     !SET .LE.0 FOR LOG MESH
      IWT=0                              !INITIAL STAT. WEIGHT FOR XSCNS
C
      IF(BFORM)READ(5,GAUSS)
C
      IF(IWT.LT.0)THEN
        CONX=1                                              ! (PI*A0**2)
        IWT=-IWT
      ENDIF
C
      IF(ILOG.LT.0.AND.E0.LE.0)E0=UNITS*TEN**ILOG     !INITIAL E FOR LOG
      IF(NGAUSS.LE.0)THEN
        IF(ILOG.GE.0)THEN
          NGAUSS=201
        ELSE
          NGAUSS=-ILOG*200
        ENDIF
      ENDIF
C
      IF(UNITS.NE.ONE)THEN                          !RYDBERGS INTERNALLY
        E0=E0/UNITS
        EMIN=EMIN/UNITS
        EMAX=EMAX/UNITS
        EWIDTH=EWIDTH/UNITS
        TPAR=TPAR/UNITS
      ENDIF
C
      IF(EWIDTH.LT.0)EWIDTH=1                           !NO MAXWELL HERE
      IF(EWIDTH.EQ.0.OR.TPAR.GT.0)THEN                  !E-DEPEND GAUSSN
        IF(TPAR.LE.0)STOP 'EWIDTH.EQ.0 REQUIRES TPAR.GT.0'
        A0=2*SQRT(TPAR)
        A0=A0*CONG
        E0=MAX(E0,TPAR)
        EWIDTH=0
        PIEA=0                              !CANNOT CHECK MESH RANGE YET
      ELSE
        A=CONG/EWIDTH
        PIEA=PIE/A
      ENDIF
C
      INQUIRE (FILE='xout',EXIST=BFORM)
      IF(BFORM)THEN
        OPEN(7,FILE='xout',STATUS='OLD',FORM='FORMATTED')    !INPUT FILE
      ELSE
        STOP 'FILE xout NOT FOUND!'
      ENDIF
C
      OPEN(8,FILE='xoutg',STATUS='UNKNOWN')           !OUTPUT ENERGY AND
C                                 GAUSSIAN CONVOLUTED COLLISION STRENGTH
      READ(7,100)IS,JS,DELE
      WRITE(8,101)IS,JS,DELE
C
      DO I=1,NMXE
        READ(7,200,END=900)E1,E2,E3,OMEGA
        EVEC(I)=E2
        VEC(I)=OMEGA
      ENDDO
      I=NMXE+1
C
 900  MXE=I-1
      IF(E1.GE.EINF)THEN                   !REMOVE INFINITE ENERGY POINT
        DEL0=0                                       !**** TEMP FIX ****
        MXE=MXE-1
      ELSE
        DEL0=E1-E2
      ENDIF
C
      IF(IWT.GT.0)THEN                         !CONVERT TO CROSS SECTION
        DO I=1,MXE
          VEC(I)=CONX*VEC(I)/(IWT*EVEC(I))
        ENDDO
      ENDIF
C
      IF(EWIDTH.EQ.0)THEN
        IF(EMAX.GT.0)THEN
          E=EMAX
        ELSE
          E=EVEC(MXE)
        ENDIF
        PIEA=PIE*A0*SQRT(E)/CONG
      ENDIF
C
      IF(EVEC(MXE).LT.EMAX+PIEA)THEN
        IF(MXE.EQ.NMXE)THEN
          WRITE(*,*)'REDUCING EMAX, INCREASE NMXE'
        ELSE
          WRITE(*,*)'REDUCING EMAX, BASED-ON xout CONTENT'
        ENDIF
        EMAX=EVEC(MXE)-PIEA
        WRITE(*,*)'EMAX=',EMAX
      ELSEIF(EMAX.LE.0)THEN                                  !USER UNSET
        EMAX=EVEC(MXE)-PIEA
      ENDIF
C
      IF(EMIN.LE.0)EMIN=E0                                   !USER UNSET
      IF(EVEC(1).GT.EMIN-PIEA.AND.VEC(1).NE.0)THEN     !NOT AT THRESHOLD
        WRITE(*,*)'INCREASING EMIN, BASED-ON xout CONTENT'
        EMIN=EVEC(1)+PIEA
        WRITE(*,*)'EMIN=',EMIN
      ENDIF
C
      NT=NGAUSS
      T=NT-1
      IF(ILOG.GT.0)THEN                                          !LINEAR
        E0=EMIN
        DEG=EMAX
      ELSE
        E0=LOG10(EMIN)
        DEG=LOG10(EMAX)
      ENDIF
      DEG=(DEG-E0)/T
C
      DO N=1,NT
        E=E0+(N-1)*DEG
        IF(ILOG.LE.0)E=TEN**E
        IF(TPAR.GT.0)EWIDTH=A0*SQRT(E)
C
        SRR=RCONVOLG(E,EWIDTH,EVEC,VEC,MXE)
C
        E2=E
        E1=E2+DEL0
        E3=E2-DELE
        WRITE(8,200)E1*UNITS,E2*UNITS,E3*UNITS,SRR
      ENDDO
C
      WRITE(0,201)IS,JS,EWIDTH*UNITS
C
 100  FORMAT(1X,I3,2X,I4,20X,F12.6)
 101  FORMAT('#',I3,' -',I4,' TRANSITION, ENERGY=',F12.6)
 200  FORMAT(1X,4(1PE14.6))
 201  FORMAT('FILE xout/g: TRANSITION',I3,' -',I4,' CONVOLUTED WITH',
     XF12.3,' FWHM GAUSSIAN')
C
      END
C
C***********************************************************************
C
      REAL*8 FUNCTION RCONVOLG(E,EWIDTH,EVEC,VEC,MXE)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C CONVOLUTE R-MATRIX COLLISION STRENGTH/CROSS SECTION WITH GAUSSIAN
C
      PARAMETER (ONE=1)
      PARAMETER (TWO=2)
C
      DIMENSION EVEC(*),VEC(*)
C
C DETERMINE CONVOLUTION ENERGY RANGE
C
      PIE=ACOS(-ONE)
      SQPIE=SQRT(PIE)
      CONG=2*SQRT(LOG(TWO))
C
      A=CONG/EWIDTH
      EMIN=E-PIE/A
      EMAX=E+PIE/A
C
C CARRY-OUT TRAPEZOIDAL RULE INTEGRATION.
C DOES **NOT** NEED CONSTANT SPACING IN ENERGY.
C CURRENTLY ASSUMES A FINE ENOUGH ENERGY MESH INPUT.
C TBD, INTERPOLATE IF A COARSER ENERGY MESH INPUT. 
C
      ANS=0
      V2=0
      DO 1 M=2,MXE
        IF(EVEC(M).LT.EMIN)GO TO 1
        H=EVEC(M)-EVEC(M-1)
        V1=V2
        T2=(EVEC(M)-E)*A
        T2=T2*T2
        V2=VEC(M)*EXP(-T2)
        ANS=ANS+H*(V2+V1)
        IF(EVEC(M).GT.EMAX)GO TO 2
   1  ENDDO
C
   2  ANS=ANS/2
C
      RCONVOLG=ANS*A/SQPIE
C
      RETURN
      END