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