Title :ELSKKT Keywords :EELS,CTEM,OPTICAL PROPERTIES Computer :11/730-785 Operating System :Vax/Vms Programing Language :Fortran IV Hardware Requirements :None Author(s) :Alex H. Buxbaum Correspondence Address :Department of Materials Engineering, Technion- Israel Institute of Technology Haifa 32000, Israel. Abstract: ELSKKT is program written for the purpose of obtaining the loss function and dielectric constant (real and imaginary parts) of materials, from Electron Energy Loss (EELS) spectra obtained in the transmission electron microscope. Usefull data generally begines around 8eV, for an instrument with a resolution of 2 eV. After correcting for multiple scattering, data is treated by using the Kramers-Kronig Analysis. Title :ELSKKT Keywords :EELS,CTEM,OPTICAL PROPERTIES Computer :11/730-785 Operating System :VAX/VMS Programing Language :Fortran IV Hardware Requirements :None Author(s) :Alex H. Buxbaum Correspondence Address :Department of Materials Engineering, Technion-Israel Institute of Technology Haifa 32000, Israel. Description: ELSKKT program was written to extract the loss function and the dielectric constant of a material from the low loss region of an EELS spectrum. This is accomplished by first deconvolving the multiply scattered spectrum (by the Fourier Log method) to obtain a single scattering distribution (SSD) from the EELS spectrum and then performing a Kramers Kronig analysis. References: Buxbaum, A.H. (1988) METHODOLOGY FOR THE STUDY OF OPTICAL PROPERTIES OF SOLIDS BY EELS WITH AN APPLICATION TO YBa2Cu3O7,Master's Thesis, Northwestern University, Evanston, Illinois. EGERTON, R., F. (1986) ELECTRON ENERGY LOSS SPECTROSCOPY IN THE ELECTRON MICROSCOPE, Plenum Press, New York The following procedure is used to compile and link ELSKKT for use on a VAX $FORTRAN ELSKKT <----------- This creates the file ELSKKT.OBJ $LINK ELSKKT <----------- This creates the file ELSKKT.EXE $RUN ELSKKT <----------- This runs the program Comments: The program reads its data from raw data file formatted to the standard EMMPDL data format. A sample data file (AL.DAT) is given at the end of the code. This data file was copied from the text by Egerton (mentioned above), to facilitate an easy comparison with the code given in Egerton's book. Example of program output using the EELS spectrum AL.DAT (included at the eind of the file ELSKKT.SRC): vmsa> RUN ELSKKT ** *** *** *** ELSKKT *** *** *** *** PLEASE MAKE SURE YOUR CAPS LOCK IS TURNED ON ENTER DATA FILE NAME (FILENAME.DAT): AL 5-OCT-90 at 21:17:16 Input Data File = AL Sample: Aluminum Spectrum = Example from Egerton Text THIS SPECTRAL FILE CONTAINS 79 CHANNELS OF DATA EXPERIMENTAL CALIBRATION CONSTANTS: OFFSET = 0.000 eV EV/CHANNEL = 1.000 ACCELERATING VOLTAGE = 100.0 kV INCID. BEAM DIVERGENCE = 5.300 mR SCATTERING ANGLE = 20.000 mR LIVE TIME = HIT TO CONTINUE IS MATERIAL A METAL (DOES IT HAVE A HIGH REFRACTIVE INDEX) [Y] ? Y DO YOU WISH TO VIEW A DETAILED OUTPUT OF ALL VALUES [N] ? Y ENTER NUMBER OF LINES OF THE RAW SPECTRUM THAT YOU WISH TO VIEW: 10 RAW SPECTRUM: 2 4 7 25 166 489 702 539 261 119 61 38 31 31 31 32 33 37 50 103 236 420 450 329 196 129 94 77 70 67 64 60 59 59 67 92 144 179 178 137 103 80 68 61 57 55 51 49 47 45 46 52 61 68 65 58 49 42 38 36 34 32 31 29 28 27 27 26 28 28 HIGHEST NUBER OF COUNTS IN ZERO LOSS PEAK = 702 NUMBER OF CHANNLES ALLOCATED TO SPECTRAL ARRAY = 512 ***THE SPECTRUM HAS BEEN EDITED FOR FURTHER PROCESSING*** HOW MANY LINES OF THE EDITED SPECTRUM DO YOU WISH TO VIEW? 10 2 4 7 25 166 489 702 539 261 119 61 38 31 31 31 32 33 37 50 103 236 420 450 329 196 129 94 77 70 67 64 60 59 59 67 92 144 179 178 137 103 80 68 61 57 55 51 49 47 45 46 52 61 68 65 58 49 42 38 36 34 32 31 29 28 27 27 26 28 28 SUBTRACT INSTRUMENT BACKGROUND?, IB= 2 COUNTS [Y] Y CHANNEL WITH MAXIMUM DATA COUNTS (ZERO LOSS PEAK) = 7 MAXIMUM DATA COUNTS (ZERO LOSS PEAK) = 702 CHANNEL OF FIRST MINIMUM AFTER THE ZERO LOSS PEAK = 14 COUNTS AT THE FIRST MINIMUM AFTER THE ZERO LOSS PEAK = 31 COUNTS IN THE LAST DATA CHANNEL OF SPECTRUM 18 HIT TO CONTINUE A GAIN CHANGE OF 1.0 WAS USED IN THIS SPECTRUM THE END OF THE SPECTRUM WILL BE FITTED TO A FUNCTION OF THE FORM A*EXP(-R), WHERE VALUES OF A AND R ARE FITTED ACCORDING TO THE DATA. R SHOULD BE -5 < R < -2. RESULTS FROM FITTING THE SPECTRUM GAVE A VALUE OF R= -4.610 CHANNEL OF GAIN CHANGE = 155 COUNTS AT GAIN CHANGE = 16.00 HOW MANY LINES OF THE UNPROCESSED SPECTRUM DO YOU WISH TO VIEW 10 UNPROCESSED SPECTRUM PRIOR TO DECONVOLUTION: 0.0 2.0 5.0 23.0 164.0 487.0 700.0 537.0 259.0 117.0 59.0 36.0 29.0 29.0 29.0 30.0 31.0 35.0 48.0 101.0 234.0 418.0 448.0 327.0 194.0 127.0 92.0 75.0 68.0 65.0 62.0 58.0 57.0 57.0 65.0 90.0 142.0 177.0 176.0 135.0 101.0 78.0 66.0 59.0 55.0 53.0 49.0 47.0 45.0 43.0 44.0 50.0 59.0 66.0 63.0 56.0 47.0 40.0 36.0 34.0 32.0 30.0 29.0 27.0 26.0 25.0 25.0 24.0 26.0 26.0 HIT TO CONTINUE ZERO LOSS PEAK STATISTICS: LOWER AND UPPER BOUND CHANNELS USED FOR FITTING THE ZERO LOSS PEAK 5 10 FWHM OF FITTED ZERO LOSS PEAK = 2.62 EV PEAK OF ZERO LOSS PEAK IS AT 7.00 EV FULL WIDTH AT TENTH MAXIMUM OF THE ZERO LOSS PEAK: 10.00 EV ENERGY WHERE THE FIRST MINIMUM AFTER THE ZERO LOSS LOSS PEAK OCCURES 14.00 EV THE ENERGY AT WHICH THE EXTRAPOLATION OF THE ZERO LOSS PEAK WILL START WILL BE 12.00, IS THIS REASONABLE? [Y] Y ZERO LOSS PEAK WILL BE EXTRAPOLATED FROM 12.00EV BY AN EXPONENTIAL CURVE AREA UNDER ZERO LOSS PEAK: 2406.3 ZERO LOSS PEAK BEFORE DECONVOLUTION: 0.00 2.00 5.00 23.00 164.00 487.00 700.00 537.00 259.00 117.00 59.00 36.00 11.69 3.79 1.23 0.40 0.13 0.04 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 HOW MANY LINES OF DECONVOLUTED SPECTRUM DO YOU WISH TO VIEW? 10 DECONVOLUTED SPECTRUM (SINGLE SCATTERING DISTRUBUTION): 0.0 0.0 -0.1 0.1 1.4 6.9 17.1 26.0 29.8 31.2 33.7 43.6 76.8 155.0 281.2 396.4 411.1 315.7 191.4 107.1 67.9 52.2 45.2 40.9 37.0 32.8 29.0 25.0 21.2 19.0 18.6 19.6 21.0 20.1 16.5 12.8 10.5 9.3 8.8 8.2 7.6 7.2 6.5 5.7 4.7 4.0 4.2 5.2 6.0 6.4 5.9 4.8 4.2 4.0 3.7 3.4 3.2 3.1 3.1 3.4 3.3 2.9 2.5 2.0 1.2 0.6 0.6 1.2 1.5 1.4 TOTAL AREA UNDER SPECTRUM = 7624.0 FOIL THICKNESS = 105.7 NM NORMALIZATION FACTOR,RK, = 9.96 UP TO WHAT ENERGY DO YOU WISH TO VIEW VALUES FOR THE DIELECTRIC FUNCTION? 20 CREATE AN OUTPUT FILE [N]? Y ENTER NAME OF OUTPUT FILE: ALOUT EV RE(1/EPS) IM(-1/EPS) EPS1 EPS2 1.000 -0.002 0.000 -453.417 -1.273 2.000 -0.021 -0.001 -47.725 -1.385 3.000 -0.055 0.001 -18.330 0.290 4.000 -0.111 0.010 -8.912 0.822 5.000 -0.189 0.052 -4.904 1.355 6.000 -0.247 0.134 -3.130 1.690 7.000 -0.271 0.208 -2.323 1.786 8.000 -0.297 0.244 -2.009 1.651 9.000 -0.375 0.260 -1.800 1.249 10.000 -0.527 0.286 -1.465 0.795 11.000 -0.798 0.377 -1.025 0.484 12.000 -1.194 0.673 -0.636 0.358 13.000 -1.559 1.379 -0.360 0.318 14.000 -1.457 2.535 -0.170 0.297 15.000 -0.446 3.618 -0.034 0.272 16.000 1.158 3.796 0.074 0.241 17.000 2.396 2.948 0.166 0.204 18.000 2.744 1.806 0.254 0.167 19.000 2.505 1.021 0.342 0.139 20.000 2.160 0.654 0.424 0.128 FORTRAN STOP Title :ELSKKT Keywords :EELS,CTEM,OPTICAL PROPERTIES Computer :11/730-785 Operating System :Vax/Vms Programing Language :Fortran IV Hardware Requirements :None Author(s) :Alex H. Buxbaum Correspondence Address :Department of Materials Engineering, Technion- Israel Institute of Technology Haifa 32000, Israel. C PROGRAM ELSKKT C C WRITTEN BY: A. BUXBAUM C C PURPOSE: OBTAIN LOSS FUNCTION AND DIELECTRIC FUNCTION C FROM AN ELECTRON ENERGY LOSS SPECTRUM. C METHOD: DECONVOLVE BY EXTRAPOLATED ZERO-LOSS PEAK C ARRAY, USING THE FOURIER-LOG METHOD. NORMALIZE C LOSS FUNCTION BY REFRACTIVE INDEX. OBTAIN DIELECTRIC C FUNCTION BY KRAMERS-KRONIG ANALYSIS. C C C REFERENCES: C C BUXBAUM, A.H. (1988), METHODOLOGY FOR THE STUDY OF C OPTICAL PROPERITES OF SOLIDS BY EELS WITH AN C APPLICATION TO YBa2Cu3O7, MASTER'S THESIS, NORTHWESTERN C UNIVERSITY, EVANSTON, ILLINOISE. C C EGERTON, R., F. (1986) ELECTRON ENERGY LOSS SPECTROSCOPY C IN THE ELECTRON MICROSCOPE. PLENUM PRESS, NEW YORK. C C DEFINE ARRAYS USED IN PROGRAM: C EC INDICATES EVERY CHANNEL ARRAY, I.E. A ONE DIMENSIONAL C ARRAY WITH CONSECUTIVE VALUES. C EOC INDICATES EVERY OTHER CHANNEL ARRAY, I.E. A ONE C DIMENSIONAL ARRAY WITH EVERY OTHER VALUE BEING ZERO C (E.G. 3,0,6,0,9,0,7....) C C ID - ORIGINAL DATA. SIZE = ND ,UP TO 1024. EC C Z - ORIGINAL DATA UP TO FIRST MINIMUM. SIZE = 100. EC C ZZ - EXTRAPOLATED IZ ARRAY WITH BACKGROUND SUBTRACTED. C SIZE = MM. EOC C D - EXTRAPOLATED SPECTRUM. SIZE = MM. EOC C BY - ZERO LOSS PEAK RANGING FROM LOWER TO UPPER BOUNDS C (IBL, IBU). TO BE SENT TO ZLPKFIT FOR A GAUSSIAN FIT. C G - GAUSSIAN CREATED FROM CALCULATIONS IN ZLPKFIT. SIZE = NN. C EC. TO BE REDEFINED. C GG - G ARRAY SHIFTED SO PEAK IS AT SAME CHANNEL AS THE MAXIMUM C DATA POINT IN THE EXPERIMENTAL ZERO LOSS PEAK. SIZE = NN.EC C G - REDEFINED AS GG. SIZE = MM. EOC C UG - ARRAY G CONVERTED TO UNIT AREA. SIZE = MM. EOC C c INTEGER*2 n INTEGER E REAL LABEL(20) CHARACTER*80 INAME COMMON /ID/ ID COMMON /RPDL/XS COMMON /D/ D COMMON /SSD/ SSD COMMON /UG/ UG,GG COMMON /ZZ/ ZZ COMMON /KKA/ PSA1,PSA2,REP,EV,AIME COMMON /ZLP/ ZLP1 C DIMENSION ID(4096) DIMENSION XS(4110) DIMENSION D(8192) DIMENSION UG(8192),GG(8192) DIMENSION ZZ(8192) DIMENSION PSA1(4096),PSA2(4096),REP(4096),EV(4096),AIME(4096) DIMENSION ZLP1(4096) c C WRITE(6,3) 3 FORMAT(/,'*** *** *** *** ELSKKT *** *** *** *** ',//, @ ' PLEASE MAKE SURE YOUR CAPS LOCK IS TURNED ON',//) C C C C *************************************************************** C READ IN EELS SPECTRUM C *************************************************************** C C XS(I) REAL ARRAY THAT CONTAINS RAW SPECTRUM C ID(I) INTEGER ARRAY THAT CONTAINS RAW SPECTRUM C NN NUMBER OF CHANNELS USED TO STORE SPECTRAL DATA AND C USED IN FOURIER ANALYSIS C ND= NUMBER DATA POINTS TO READ C EPC= ENERGY PER CHANNEL C E0= ACCELERATING VOLTAGE C BETA= SPECTROMETER CONVERGENCE ANGLE C RI= REFRACTIVE INDEX OF SPECIMEN C 1 CONTINUE C DO 2 I=1,4096 XS(I)=0. 2 ID(I)=0 C CALL RPDL(ND,EPC,E0,BETA,IOUTNAME) C C TRANSFER DATA FROM ARRAY XS (REAL VALUES) INTO ARRY ID C (ITNEGER VALUES) DO 61 I=1,ND-1 61 ID(I)=IFIX(XS(I+1)) C C SUPPLY A DEFAULT VALUE FOR THE REFRACTIVE INDEX RI = 1000 C WRITE(6,55) 55 FORMAT(' IS MATERIAL A METAL (DOES IT HAVE A HIGH REFRACTIVE @ INDEX) [Y] ?') READ(6,67)ANSWER4 IF(ANSWER4.EQ.'N') THEN WRITE(6,56) 56 FORMAT(' ENTER REFRACTIVE INDEX...') READ(6,*)RI ELSE CONTINUE ENDIF C GIVE THE OPTION OF DUMPING ALL CALCULATED VALUES IN THIS PRPROGRAM C TO THE TERMINAL IPRINT=-1 64 WRITE(6,66) 66 FORMAT(' DO YOU WISH TO VIEW A DETAILED OUTPUT OF ALL @ VALUES [N] ? ') READ (6,67) ANSWER1 67 FORMAT(1A1) IF(ANSWER1.EQ.'Y') IPRINT=+1 IF(IPRINT.LT.0) GOTO 99 WRITE(6,68) READ (6,69) ILAST 68 FORMAT(/,' ENTER NUMBER OF LINES OF THE RAW SPECTRUM THAT YOU @ WISH TO VIEW: ') 69 FORMAT(1I3) WRITE(6,70) 70 FORMAT(' RAW SPECTRUM:') WRITE(6,75)(ID(I),ID(I+1),ID(I+2),ID(I+3),ID(I+4),ID(I+5), @ ID(I+6),I=1,ILAST*7,7) 75 FORMAT(7I8) 99 CONTINUE C C ************************************************************** C FIND FIRST AND LAST DATA POINTS IN RAW SPECTRUM. C GET RID OF CHANNELS WITH NO COUNTS ON GAIN SIDE C OF ZERO LOSS PEAK. C CALL RANGE(ND,NN) WRITE(6,100)NN 100 FORMAT('0NUMBER OF CHANNLES ALLOCATED TO SPECTRAL ARRAY = @ ',I5) WRITE(6,125) 125 FORMAT('0***THE SPECTRUM HAS BEEN EDITED FOR FURTHER @PROCESSING***') IF(IPRINT.LT.0) GOTO 199 WRITE(6,126) 126 FORMAT(/,'0HOW MANY LINES OF THE EDITED SPECTRUM @ DO YOU WISH TO VIEW?') READ(6,127)ILAST 127 FORMAT(1I3) DO 130 I=1,ILAST*7,7 130 WRITE(6,140)ID(I),ID(I+1),ID(I+2),ID(I+3),ID(I+4),ID(I+5) @ ,ID(I+6) 140 FORMAT(1X,7I8) 199 CONTINUE C ************************************************************** C FIND FIRST MAXIMUM AND MINIMUM IN DATA CALL SHAPE(NN,ND,IB,NOISE,NZ,MMIN ,MFIN,MSEP,MM) IF(IPRINT.LT.0) GOTO 299 WRITE(6,232)NZ 232 FORMAT('0CHANNEL WITH MAXIMUM DATA COUNTS (ZERO LOSS PEAK) = ' @ ,I3) WRITE(6,234)ID(NZ) 234 FORMAT('0MAXIMUM DATA COUNTS (ZERO LOSS PEAK) = ',I8) WRITE(6,237)MMIN 237 FORMAT('0CHANNEL OF FIRST MINIMUM AFTER THE ZERO LOSS PEAK = @ ',I3) WRITE(6,239)ID(MMIN) 239 FORMAT('0COUNTS AT THE FIRST MINIMUM AFTER THE ZERO LOSS PEAK = @ ',I6) WRITE(6,242)ID(ND) 242 FORMAT ('0COUNTS IN THE LAST DATA CHANNEL OF SPECTRUM',I5) 299 CONTINUE C ************************************************************** WRITE(6,300) 300 FORMAT(' HIT TO CONTINUE ') READ(6,301)X 301 FORMAT(80A1) C ELIMINATE GAIN CHANGE CALL GAIN(ND,NOISE,IB,MMIN,GC,MG) IF(IPRINT.LT.0) GOTO 321 WRITE(6,315)GC 315 FORMAT('0A GAIN CHANGE OF ',F5.1,' WAS USED IN THIS SPECTRUM') 321 CONTINUE C ************************************************************** C GENERATE SPECTRUM FOR ANALYSIS C EXRAPOLATE LAST DATA CHANNEL TO ZERO C STATEMENTS 400-499 CALL SPEC(ND,MG,MFIN,IB,MM,GC,TAREA) IF(IPRINT.LT.0) GOTO 499 WRITE(6,460)MG 460 FORMAT('0CHANNEL OF GAIN CHANGE = ',I5) WRITE(6,455)D(MG) 455 FORMAT( '0COUNTS AT GAIN CHANGE = ',F10.2) WRITE(6,456) 456 FORMAT(/,' HOW MANY LINES OF THE UNPROCESSED SPECTRUM @ DO YOU WISH TO VIEW') READ(6,457)ILAST 457 FORMAT(1I3) WRITE (6,465) 465 FORMAT ('0UNPROCESSED SPECTRUM PRIOR TO DECONVOLUTION:') DO 470 I=1,ILAST*14,14 470 WRITE (6,475)D(I),D(I+2),D(I+4),D(I+6),D(I+8),D(I+10),D(I+12) 475 FORMAT(1X,7F10.1) 499 CONTINUE C ************************************************************** C FIT GAUSSIAN TO ZERO LOSS PEAK BY LEAST SQUARES USING C THE SUBROUTINE ZLPKFIT. C CALL ZEROFIT(MMIN,NN,NZ,AMP,SIGFIT,IFWTM,IBL) WRITE(6,481)(2.*SIGFIT*0.8326)*EPC 481 FORMAT('0FWHM OF FITTED ZERO LOSS PEAK = ',F10.2,' EV') C C ************************************************************* C C ************************************************************* C SET UP ZERO LOSS SPECTRUM C A0 IS COUNTS UNDER ZERO LOSS PEAK C ZZ IS ZERO LOSS SPECTRUM (EVERY OTHER CHANNEL ARRAY) C C CALL ZEROEXT(SIGFIT,IBL,EPC,IFWTM,NZ,MMIN,IB,MSEP,MM,NN,A0) C STATEMENTS 1200-1299 C DEFINES ZERO LOSS ARRAY FROM INELASTIC SPECTRUM BY C EXTRAPOLATING ZERO LOSS TO ZERO FROM A SPECIFIED CHANNEL C (COMMONLY FROM THE CHANNEL WHERE THE ZERO LOSS PEAK IS C 10% OF ITS MAXIMUM). C OUTPUT IN ARRAY ZZ (EVERY OTHER CHANNEL ARRAY) C C ************************************************************** C SET UP A UNIT AREA GAUSSIAN WITH SIGMA OF EXPERIMENTAL ZERO LOSS PEAK C GAUSSIAN IN ARRAY GG(J) (EVERY OTHER CHANNEL ARRAY) C UNIT AREA GAUSSIAN IN ARRAY UG(J) (EVERY OTHER CHANNEL ARRAY) C C EQUATE THE STANDARD DEVIATION OF THE EXPERIMENTAL ZERO LOSS C PEAK (SIGFIT) TO THE STANDARD DEVIATION OF THE GAUSIAN TO BE C USED FOR RECONVOLUTION (SIG) C SIG=SIGFIT C CALL GAUS(NN,NZ,MMIN,SIG,AMP) C ************************************************************* C DECONVOLVE SPECTRUM BY FOURIER-LOG METHOD C DATA IN ARRAY D (EVERY OTHER CHANNEL ARRAY) C ZERO LOSS PEAK IN ARRAY ZZ (EVERY OTHER CHANNEL ARRAY) C UNIT AREA GAUSSIAN IN ARRAY UG (EVERY OTHER CHANNEL ARRAY) C C DECONVOLVED SPECTRUM IS WRITTEN BACK IN ARRAY D CALL DECON(NN,MM,A0) IF(IPRINT.LT.0) GOTO 799 WRITE (6,705) 705 FORMAT ('0HOW MANY LINES OF DECONVOLUTED SPECTRUM DO YOU WISH @ TO VIEW?') READ(6,*)DLAST WRITE (6,710) 710 FORMAT ('0DECONVOLUTED SPECTRUM (SINGLE SCATTERING @DISTRUBUTION):') DO 720 J=1,DLAST*14,14 720 WRITE (6,725)D(J),D(J+2),D(J+4),D(J+6),D(J+8),D(J+10),D(J+12) 725 FORMAT(1X,7(F10.1)) 799 CONTINUE C *********************************************************** C KRAMERS KRONIG ANALYSIS C NORMALIZE ARRAY D TO OBTAIN THE LOSS FUNCTION C PERFORM KRAMERS KRONIG ANALYSIS TO OBTAIN Re(1/E) C CALCULATE COMPLEX DIELECTRIC FUNCTION CALL KKA(TAREA,MM,NN,NZ,MSEP,EPC,E0,BETA,RI,A0,ANISO) C WRITE(6,805) 805 FORMAT(/,' UP TO WHAT ENERGY DO YOU WISH TO VIEW VALUES FOR @THE DIELECTRIC FUNCTION?') READ(6,*)DLAST WRITE(6,62) 62 FORMAT(' CREATE AN OUTPUT FILE [N]? ') READ (6,67) ANSWER2 IF(ANSWER2.EQ.'Y') THEN WRITE(6,58) 58 FORMAT('0ENTER NAME OF OUTPUT FILE:') READ(6,59)INAME 59 FORMAT(A80) ENDIF C IF(ANSWER2.EQ.'Y') THEN C C OPEN (UNIT=12,NAME=INAME,TYPE='NEW',ACCESS='SEQUENTIAL') DO 897 I=1,IFIX(DLAST/EPC) WRITE(12,899) EV(I),REP(I),AIME(I+1),PSA1(I),PSA2(I) 897 CONTINUE CLOSE(UNIT=12) ENDIF WRITE (6,898) 898 FORMAT(/' ',5X,'EV',4X,'RE(1/EPS)',2X,'IM(-1/EPS)',4X,'EPS1', 1 7X,'EPS2'/) WRITE (6,900)(EV(I),REP(I),AIME(I+1),PSA1(I),PSA2(I) @ ,I=1,IFIX(DLAST/EPC)) 899 FORMAT(F8.3,2X,4(1X,F11.4),1X,E9.1) 900 FORMAT(5F11.3) C C C****************************************************************** C END OF MAIN PROGRAM C STOP C END C C******************************************************************** C SUBROUTINES: C******************************************************************** C RANGE (ND,NN) C C SUBROUTINE RANGE(ND,NN) C C THIS SUBROUTINE DEFINES THE USEFUL PART OF A RAW SPECTRUM C AND STORES IT IN ARRAY ID. THE DATA CONSIDERED 'USEFUL' C IS FROM THE LAST MINIMUM ON THE GAIN SIDE OF THE ZERO LOSS PEAK C (WHERE THE ZERO LOSS PEAK HAS A POSITIVE SLOPE ALL THE WAY C TO ITS MAXIMUM) UNTIL THE LAST DATA POINT COLLECTED. C COMMON /ID/ ID DIMENSION ID(4096) C SET FIRST POINT IN THE RAW DATA ARRAY ID TO LAST MINIMUM ON C THE GAIN SIDE OF THE ZERO LOSS PEAK (IN CHANNEL IDIP) OR C SET FIRST POINT IN ARRAY ID TO FIRST NON-ZERO DATA POINT C C LOCATE GAIN CHANGE (CHANNEL IGAIN), THEN DEFINE PEAK OF ZERO C LOSS PEAK AS THE CHANNEL WITH MOST COUNTS PRIOR TO GAIN CHANGE C (CHANNEL IPEAK). C ID - RAW DATA ARRAY (INTEGER, SINGLE CHANNELED C ND - CHANNEL OF LAST DATA POINT C IGAIN - CHANNEL OF GAIN CHANGE C C DO 110 I3=10,ND 110 IF (ID(I3).GT.10*ID(I3-1).AND.ID(I3)*ID(I3-1).NE.0)GOTO 120 120 IGAIN=I3 IPEAK1=1 DO 130 I=1,IGAIN-1 130 IF(ID(I).GT.ID(IPEAK1)) IPEAK1=I WRITE(6,135)ID(IPEAK1) 135 FORMAT('0HIGHEST NUBER OF COUNTS IN ZERO LOSS PEAK = @',I10) ICNT=0 C DEFINE FIRST CHANNEL OF USEFUL DATA TO BE THE FIRST MINIMUM ON C THE GAIN SIDE OF THE ZERO LOSS PEAK. C C C FOR A MINIMUM OF A SINGLE CHANNEL DO 140 I=1,IPEAK1 IF(ID(I).LT.ID(I-1).AND.ID(I).LT.100) IDIP=I 140 CONTINUE C IF(IDIP.NE.0) THEN ICNT=IDIP GO TO 150 ELSE ICNT=0 DO 145 I=1,IPEAK1 ICNT=ICNT+1 145 IF(ID(I).NE.0) GO TO 150 ENDIF 150 CONTINUE C C SHIFT SPECTRUM SO FIRST USEFUL DATA IS IN THE FIRST CHANNEL DO 160 I=1,ND+ICNT 160 ID(I)=ID(I+ICNT-1) ND=ND+1-ICNT C C ND IS LAST RAW DATA POINT C SET ND TO LAST NON ZERO VALUE DO 170 I=1,ND 170 IF(ID(I).EQ.0) GO TO 175 175 ND=I-1 C C SET SIZE OF SPECTRAL ARRAY C SIZE MUST ME A MULTIPLE OF 2**I (DUE TO THE NATURE OF THE C FOURIER TRANSFORM ROUTINE). DO 180 I=9,12 IF(ND+50.LE.2**I) GO TO 190 180 CONTINUE 190 NN=2**I IF(NN.GT.4096) NN=4096 RETURN END C ************************************************************ C C THIS ROUTINE EXAMINES THE INPUT SPECTRAL ARRAY (ID) C AND DOES THE FOLLOWING: C 1) SUBTRACTS BACKGROUND, IB (WHICH IS FIRST DATA IN ARRAY) C 2) LOCATES PEAK OF ZERO LOSS, CHANNEL NZ C 3) FINDS FIRST MINIMUM AFTER ZERO LOSS, CHANNEL MMIN C 4) DEFINES VALUES USED IN EVERY OTHER CHANNEL ARRAY (EOCA): C MSEP IS FIRST MINIMUM AFTER ZERO LOSS PEAK C MFIN IS LAST COLLECTED DATA POINT C MM IS SIZE OF SPECTRAL ARRAY C SUBROUTINE SHAPE(NN,ND,IB,NOISE,NZ,MMIN,MFIN,MSEP,MM) COMMON /ID/ ID DIMENSION ID(4096) C C CONSIDER FIRST CHANNEL TO BE BACKGROUND NOISE. THIS WILL C CAUSE PROBLEMS IF SPECTRUM HAS LOW SIGNAL. IB=ID(1) WRITE(6,200)ID(1) 200 FORMAT('0SUBTRACT INSTRUMENT BACKGROUND?, IB=',I4,' COUNTS [Y]') READ(6,201)ANS 201 FORMAT(1A1) IF(ANS.EQ.'N') IB=0 C FIND FIRST MAXIMUM IN DATA ARRAY ID(I): C NZ IS CHANNEL OF MAXIMUM DATA POINT; ID(NZ) IS MAX DATA C IF DATA IS LESS THAN NOISE, SKIP IT AND MOVE ON NOISE=ID(1)+ID(2)+IFIX(SQRT((FLOAT(ID(1)+ID(2))))) DO 202 I1=3,ND IF (ID(I1).LT.NOISE) GO TO 202 C IF (ID(I1).LT.ID(I1-1).AND.ID(I1-1).LT.ID(I1-2))GO TO 205 IF(ID(I1).GT.ID(I1-1).AND.ID(I1).GT.ID(I1+1).AND.ID(I1).GT. @ ID(I1-2).AND.ID(I1).GT.ID(I1+2).AND.ID(I1).GT.ID(I1-3).AND. @ ID(I1).GT.ID(I1+3)) GO TO 205 202 CONTINUE C 205 I1=I1-2 !CORRECT FOR EXAMINING TWO CHANNELS BEYOND PEAK 205 NZ=I1 C FIND FIRST MINIMUM AFTER ZERO LOSS PEAK, ID(I2-1) DO 210 I2=NZ+2,ND IF(ID(I2).GT.ID(I2-2).AND.ID(I2-1).GT.ID(I2-2))GO TO 215 210 CONTINUE 215 I2=I2-2 !CORRECT FOR EXAMINING TWO CHANNELS BEYOND MIN C MMIN IS FIRST MINIMUM AFTER ZERO LOSS PEAK MMIN=I2 C SET MMIN TO BE FIRST OF TWO POINTS IF C TWO CONSECUTIVE DATA ARE EQUAL IF(ID(MMIN).EQ.ID(MMIN-1)) MMIN=MMIN-1 C DEFINE VALUES USED IN EVERY OTHER CHANNEL ARRAY (EOCA) C C MSEP IS FIRST MINIMUM AFTER ZERO LOSS PEAK MSEP=2*MMIN-1 C MFIN IS LAST COLLECTED DATA POINT MFIN=2*ND-1 C MM IS SIZE OF SPECTRAL ARRAY MM=2*NN C RETURN END C C *********************************************************** C THIS ROUTINE ELIMINATES GAIN CHANGE IN SPECTRUM C IT EXAMINES THE SPECTRAL ARRAY (ID) AND DOES THE FOLLOWING: C 1) FINDS THE CHANNEL OF GAIN CHANGE (WHERE TWO CONSECUTIVE C DATA POINTS DIFFER BY AN ORDER OF MAGNITUDE). C THAT CHANNEL IS DEFINED AS 'MG' . C 2) SETS CHANNEL ON LOW END OF GAIN CHANGE TO TWICE C THE VALUE OF ONE CHANNEL BEFORE IT LESS THE VALUE OF C TWO CHANNEL BEFORE IT. C 3) SETS VALUE OF GAIN CHANGE FACTOR (GC) TO THE RATIO OF THE C VALUE AT HIGH END OF GAIN CHANGE TO TWICE THE VALUE AT LOW C END OF GAIN CHANGE LESS THE VALUE PREVIOUS TO THAT LESS C THE BACKGROUND, IB. C C SUBROUTINE GAIN(ND,NOISE,IB,MMIN,GC,MG) COMMON /ID/ ID DIMENSION ID(4096) C FIND GAIN CHANGE, IF ANY DO 300 I3=MMIN,ND IF (ID(I3).LE.NOISE) GO TO 300 IF (ID(I3).LT.10*ID(I3-1)) GO TO 300 ID(I3-1)=2*ID(I3-2)-ID(I3-3) ID(I3)=2*ID(I3+1)-ID(I3+2) GC=FLOAT(ID(I3))/FLOAT(2*ID(I3-1)-ID(I3-2)-IB) GO TO 310 300 GC=1. I3=ND 310 CONTINUE C MG IS POINT OF GAIN CHANGE MG=2*I3-1 C RETURN END C C ********************************************************** C THIS ROUTINE TAKES THE INPUT SPECTRAL ARRAY (ID), AND PREPARES C A SPECTRUM FOR DECONVOLUTION, IN ARRAY D, BY THE FOLOWING STEPS: C C 1) COPIES DATA FROM ARRAY ID TO ARRAY D FROM FIRST CHANNEL TO C POINT OF GAIN CHANGE, AT THE SAME STEP IT SUBTRACTS THE C BACKGROUND, IB. C 2) COPIES DATA FROM ARRAY ID FROM GAIN CHANGE TO LAST COLLECTED C DATA POINT IN SPECTRUM. AT THE SAME STEP, THE DATA IS DEVIDED C BY THE GAIN CHANGE FACTOR, GC. C 3) FITTING PARAMETERS ARE DETERMINED AT THE END OF COLLECTED DATA C FOR EXTRAPOLATION. FITTING PARAMETERS ARE BASED ON VALUES OF C LAST TEN DATA POINTS. C 4) DATA IS EXTRAPOLATED FROM LAST COLLECTED DATA POINT DOWN TO C CHANNEL NN, BY THE FORM A*E**R WHERE R IS A NEGATIVE C FITTING PARAMETER. SUBROUTINE SPEC(ND,MG,MFIN,IB,MM,GC,TAREA) COMMON /ID/ ID COMMON /D/ D DIMENSION ID(4096) DIMENSION D(8192) C C C C TRANSFER SPECTRUM INTO ARRAY D C FROM FIRST CHANNEL TO GAIN CHANGE, MG DO 400 J=1,MG,2 D(J)=FLOAT(ID((J+1)/2))-FLOAT(IB) 400 CONTINUE C TRANSFER SPECTRUM INTO D(J) FROM MG TO MFIN IF(MG.EQ.MFIN) GO TO 420 DO 420 J=MG,MFIN,2 D(J)=FLOAT(ID((J+1)/2))/GC 420 CONTINUE C SET ND TO LAST NON ZERO VALUE DO 425 J=3,MFIN,2 IF(D(J).LE.0) THEN WRITE(6,421)J,MFIN,D(J) 421 FORMAT(1X,I5,2X,I7,2X,F10.4) MFIN=J-2 ND=(MFIN+1)/2 GO TO 430 ENDIF 425 CONTINUE 430 CONTINUE C C EXTRAPOLATE THE SPECTRUM TO ZERO AT J=MFIN BY FITTING TO A C FUNCTION OF THE FORME A*E^(-R) C EXTRAPOLATE SPECTRUM BY CALLING SPECEXT WHICH RETURNS C EXTTRAPOLATION PARAMETERS A AND R CALL SPECEXT(MFIN,A,R) DO 450 I=MFIN+2,MM,2 D(I)=A*FLOAT((I+1)/2)**R 450 D(I+1)=0. D(MM)=0. C C FIND AREA UNDER SPECTRUM TAREA=0. DO 460 J=1,MM,2 460 TAREA=TAREA+D(J) RETURN END C C ******************************************************************* C SUBROUTINE ZEROEXT C DEFINES A ZERO LOSS PEAK FROM THE INELASTIC SPECTRUM C EXTRAPOLATES ZERO LOSS PEAK OF INELASTIC SPECTRUM SMOOTHLY C TO ZERO C SUBROUTINE ZEROEXT(SIGFIT,IBL,EPC,IFWTM,NZ,MMIN,IB @ ,MSEP,MM,NN,A0) COMMON/D/ D COMMON/ZZ/ ZZ COMMON/ZLP/ZLP1 DIMENSION D(8192),ZZ(8192),ZLP1(4096) C C DEFINE ZERO LOSS PEAK ZZ(E) C C COPY ZERO LOSS REGION OF ID ARRAY UPTO A CHANNEL SELECTED BY C THE USER (IX2). DATA IS EXTRAPOLATED FROM IX2 TO MM BY ONE OF C THREE EXRAPOLATION METHODS: A)EXPONENTIAL B)LINEAR OR C)GAUSSIAN. C SIG1=SIGFIT WRITE(6,516)FLOAT(NZ)*EPC 516 FORMAT('0PEAK OF ZERO LOSS PEAK IS AT',F6.2,' EV') WRITE(6,515)FLOAT(IFWTM)*EPC 515 FORMAT('0FULL WIDTH AT TENTH MAXIMUM OF THE ZERO LOSS PEAK:' @,F7.2,' EV') WRITE(6,510)FLOAT(MMIN)*EPC 510 FORMAT('0ENERGY WHERE THE FIRST MINIMUM AFTER THE ZERO LOSS @ LOSS PEAK OCCURES ',F7.2,' EV') WRITE(6,520)FLOAT(IFWTM+MMIN)*EPC/2. 520 FORMAT('0THE ENERGY AT WHICH THE EXTRAPOLATION OF THE ZERO' @,/,'0LOSS PEAK WILL START WILL BE ',F7.2,', IS THIS REASONABLE? @ [Y]') READ(6,521)ANSWER3 521 FORMAT(1A1) IF(ANSWER3.EQ.'N') THEN WRITE(6,522) 522 FORMAT('0ENTER THE ENERGY AT WHICH THE EXTRAPOLATION OF THE @ZERO',/,'0LOSS PEAK WILL START:') READ(6,*)EVS ELSE EVS=FLOAT(IFWTM+MMIN)*EPC/2. ENDIF IX1=IFIX(EVS/EPC) IX2=2*IX1-1 C C DATA IS COPIED INTO ZZ ARRAY FROM FIRST CHANNEL TO IX2 C DO 525 I=1,IX2 525 ZZ(I)=D(I) C C 530 WRITE(6,535)EVS 535 FORMAT('0ZERO LOSS PEAK WILL BE EXTRAPOLATED FROM ',F10.2 @ ,'EV BY AN EXPONENTIAL CURVE') CALL EXPEXT(IX1,IX2,NN) C C C CALCULATE AREA UNDER ZERO LOSS PEAK A0=0. DO 540 I=1,MM,2 540 A0=A0+ZZ(I) WRITE(6,545)A0 545 FORMAT('0AREA UNDER ZERO LOSS PEAK:',F10.1) WRITE(6,550) 550 FORMAT ('0ZERO LOSS PEAK BEFORE DECONVOLUTION:') DO 555 I=1,MM,14 IF(I.GT.NZ.AND.ZZ(I).EQ.0.0) GOTO 561 555 WRITE(6,560) ZZ(I),ZZ(I+2),ZZ(I+4),ZZ(I+6),ZZ(I+8), @ ZZ(I+10),ZZ(I+12) 560 FORMAT(7(F10.2)) C C 561 CONTINUE C RETURN END C C ***************************************************************** C C EXTRAPOLATE ZZ TO ZERO USING AN EXPONENTIAL CURVE SUBROUTINE EXPEXT(IX1,IX2,NN) COMMON /ZZ/ ZZ DIMENSION ZZ(8192) C DY1=(ZZ(IX2)-ZZ(IX2-2))/1. DY2=(ZZ(IX2-2)-ZZ(IX2-4))/1. DY=(DY1+DY2)/2. BB=DY/ZZ(IX2) C DO 500 I=IX1,NN ZZ(2*I-1)=ZZ(IX2)*EXP(BB*(FLOAT(I)-FLOAT(IX1))) 500 ZZ(2*I)=0. RETURN END C ******************************************************************* C C THIS ROUTINE CREATES A UNIT AREA GAUSSIAN FOR C RECONVOLUTION. SUBROUTINE GAUS(NN,NZ,MMIN,SIG,AMP) COMMON /UG/ UG,GG COMMON /ZZ/ ZZ DIMENSION G(4096), GG(8192), BX(4096),BY(4096) DIMENSION UG(8192) DIMENSION ZZ(8192) INTEGER E C GENERATE GAUSSIAN FROM E=1 TO NN CENT=FLOAT(NN/2) GSUM=0. DO 620 E=1,NN GE=FLOAT(E) G(E)=AMP*EXP(-(GE-CENT)**2/(2*SIG*SIG)) 620 GSUM=GSUM+G(E) CALL SHIFT(G,NN,NN/2) C GAUSSIAN CENTERED AT FIRST CHANNEL C TRANSFER ARRAY GG TO DOUBLE SPACED ARRAY, G DO 635 E=1,NN GG(2*E-1)=G(E) UG(2*E-1)=G(E)/(AMP*SQRT(2.*3.14159)*SIG) UG(2*E)=0. 635 GG(2*E)=0. C RETURN END C************************************************************ C C ******************************************************* C C SPECTRUM IN ARRAY D IS DECONVOLVED. C THE ARRAYS D,ZZ,AND UG ARE FORWARD FOURIER C TRANSFORMED AND DECONVOLUTION TAKES THE FORM: C d = I0 ug ALOG(d/zz) C WHERE LOWER CASE DENOTES FOURIER TRANSFORMED ARRAYS C NOTE: THE ARRAY IN REAL SPACE ARE EVERY OTHER CHANNEL C ARRAYS, WITH ZERO IN EVERY EVEN CHANNEL. AFTER FOURIER C TRANFORMATION THE COSINE AND SINE COEFFICIENTS ARE C STORED IN ODD AND EVEN CHANNELS. UPON REVERSE FOURIER C TRANSFORMATION, THE EVEN CHANNELS SHOULD ONCE AGAIN C BE ZERO (I.E. NO IMAGINARY PARTS). SUBROUTINE DECON(NN,MM,A0) COMMON /D/ D COMMON /ZZ/ ZZ COMMON /UG/ UG DIMENSION D(8192) DIMENSION UG(8192) DIMENSION ZZ(8192) complex t1,t2,t3,T4,T5 C C FORWARD FORIER TRANSFORM CALL FFT(NN,+1,D) CALL FFT(NN,+1,ZZ) CALL FFT(NN,+1,UG) C DO 700 I=1,MM,2 C t2=cmplx(d(i),d(i+1) ) t3=cmplx(zz(i),zz(i+1)) t1=clog(t2/t3) C t4=cmplx(UG(i),UG(i+1)) T5=A0*t4*t1 d(i)=REAL(T5)/FLOAT(NN) d(i+1)=AIMAG(T5)/FLOAT(NN) ZZ(I)=REAL(T3)/FLOAT(NN) ZZ(I+1)=AIMAG(T3)/FLOAT(NN) UG(I)=REAL(T4)/FLOAT(NN) UG(I+1)=AIMAG(T4)/FLOAT(NN) 700 CONTINUE C C REVERSE FOURIER TRANSFORM: CALL FFT(NN,-1,D) CALL FFT(NN,-1,ZZ) CALL FFT(NN,-1,UG) C RETURN END C *********************************************************** C KRAMERS KRONIG ANALYSIS SUBROUTINE SUBROUTINE KKA(TAREA,MM,NN,NZ,MSEP,EPC,E0,BETA,RI,A0,ANISO) COMMON /KKA/ PSA1,PSA2,REP,EV,AIME COMMON /RE/ RE COMMON /D/ D COMMON /SSD/ SSD DIMENSION PSA1(4096),PSA2(4096),REP(4096),EV(4096),AIME(4096) DIMENSION SSD(4096) DIMENSION D(8192) DIMENSION RE(4096) C C C STORE SINGLE SCATTERING DISTRIBUTION IN ARRAY SSD DO 800 J=1,MM,2 800 SSD((J+1)/2)=D(J) C C WRITE(6,850)TAREA 850 FORMAT(/,' TOTAL AREA UNDER SPECTRUM = ',F10.1) C C TGT IS "TWO GAMMA TAU", IS A RELATIVISTIC FACTOR USED C IN CALCULATING THETA-E. SEE EGERTON (1986, PG. 375). TGT=E0*(1022.12+E0)/(511.06+E0) 855 SUM=0. AREA=0. C APPLY APATURE CORRECTION APC DO 860 J=3,MM,2 Eng=EPC*FLOAT(J-1)/2. APC=ALOG(1.+(BETA*TGT/Eng)**2) D(J)=D(J)/APC 860 SUM=SUM+D(J)/Eng RK=SUM/1.571/(1.-1./RI/RI)*EPC C THICKNESS IN NANOMETERS (TNM): TNM=332.5*RK/(A0*EPC)*E0*(1.+E0/1022.)/(1.+E0/511.)**2 C "T OVER LAMDA" (TOL): tol=alog(TAREA/A0) C MEAN FREE PATH (RLAM): RLAM=TNM/TOL C APPLY NORMALIZATION FACTOR RK, STORING IM(-1/EPS) IN ARRAY AIME DO 870 J=1,MM,2 D(J)=D(J)/RK D(J+1)=0. AIME((J+1)/2)=D(J) 870 CONTINUE C C C C C KRAMERS KRONIG ANALYSIS USING FFT C REF: EGERTON (1986, PG. 365) CALL FFT(NN,+1,D) C TRANSFER SINE COEFFICIENTS TO COSINE LOCATION DO 875 J=1,MM,2 D(J)=2.*D(J+1)*FLOAT((J-NN)/IABS(J-NN))/FLOAT(NN) 875 D(J+1)=0. C CALL FFT(NN,-1,D) C DO 880 J=1,NN,2 C CORRECT FOR TAIL OF REFLECTED PART OF EVEN FUNCTION D(J)=D(J)+1.-D(NN-1)/2.*(FLOAT(NN-1)/FLOAT(MM-J))**2 880 D(NN+J)=1.+D(NN-1)*(FLOAT(NN-1)/FLOAT(NN+J))**2/2. C C C TAU (T) IS RELATIVISTIC ENERGY C REF: EGERTON (1986, PG. 375) T=E0*(1.+E0/2./511.06)/(1.+E0/511.06)**2 RKO=2590.*(1.+E0/511.06)*SQRT(2.*T/511.06) DO 885 J=2,NN C ENERGY LOSS (EV) EV(J-1)=FLOAT(J-1)*EPC C REAL PART OF 1/EPS (REP) REP(J-1)=D(2*J-1) DEN=D(2*J-1)*D(2*J-1)+AIME(J)*AIME(J) C REAL PART OF DIELECTRIC CONSTANT (EPS1) PSA1(J-1)=D(2*J-1)/DEN C IMAGINARY PART OF DIELECTRIC CONSTANT (EPS2) PSA2(J-1)=AIME(J)/DEN C THETA-E, CHARACTERISTIC SCATTERING ANGLE THETAE=EV(J)/TGT 885 CONTINUE WRITE (6,895) TNM,RK 895 FORMAT(/,' FOIL THICKNESS = ',F8.1,' NM',1X,'NORMALIZATION @ FACTOR,RK, = ',F7.2) 865 CONTINUE C RETURN END c C ********************************************************* C SUBROUTINE TO PERFORM FAST FOURIER TRANSFORMATION C REF EGERTON (1986, PP. 365 - 366) c c c SUBROUTINE FFT(NN,ISIGN,D) DIMENSION D(8192) N=2*NN J=1 DO 5 I=1,N,2 IF (I-J)1,2,2 1 TR=D(J) TI=D(J+1) D(J)=D(I) D(J+1)=D(I+1) D(I)=TR D(I+1)=TI 2 M=N/2 3 IF (J-M)5,5,4 4 J=J-M M=M/2 IF (M-2)5,3,3 5 J=J+M MM=2 6 ML=MM-N IF (ML)7,10,10 7 IS=2*MM TH=6.283185/FLOAT(ISIGN*MM) ST=SIN(TH/2) W1=-2.*ST*ST W2=SIN(TH) WR=1. WI=0. DO 9 M=1,MM,2 DO 8 I=M,N,IS J=I+MM TR=WR*D(J)-WI*D(J+1) TI=WR*D(J+1)+WI*D(J) D(J)=D(I)-TR D(J+1)=D(I+1)-TI D(I)=D(I)+TR 8 D(I+1)=D(I+1)+TI TR=WR WR=WR*W1-WI*W2+WR 9 WI=WI*W1+TR*W2+WI MM=IS GO TO 6 10 RETURN END C ******************************************************** C SUBROUTINE ZLPKFIT (BX,BY,IBL,IBU,AMP,EBAR,SIG) C SOURCE: NESTOR ZALUZEC C DIMENSION BX(4096),BY(4096) C C BX = ARRAY CONTAINING THE X AXIS DATA C BY = ARRAY CONTAINING THE Y AXIS DATA C C SUBROUTINE TO FIT A GAUSSIAN PEAK TO DATA USING A LEAST C SQUARES ANALYSIS THE LEAST SQUARES FIT ONLY ITERATES C THE PARAMETERS OF AMPLITUDE AND WIDTH THE ENERGY CENTROID C IS ASSUMED TO REMAIN CONSTANT C C DEF: GAUSSIAN(EI) = AMP * EXP (-0.5 *((EI-EBAR)/SIG)**2.) C C AMP = AMPLITUDE C EBAR = CENTROID ENERGY C SIG = STANDARD DEVIATION (WIDTH) C EI = CHANNEL (ENERGY) VARIABLE C C C THE FIT IS MADE OVER THE ROI WHICH WAS SELECTED BY USER C C C COMMON /PARAM1/ OFFSET,EVCH,A,R,YY,YN,EX(100),COE(7) COMMON /PARAM2/ IBFLG,IZERO,IBGND,IROI,NEL DOUBLE PRECISION SUM(5),X,Y C C C DO LINEAR LEAST SQUARES FIT OF ROI C CALL LSTSQR(1,,,,,,SUM) !USE BRESSLER'S ROUTINE HERE NUMPTS=IBU-IBL+1 C C SEND MODIFIED DATA TO LSTSQR ROUTINE C C C WE DON'T NEED INITIAL PARAMETERS FOR AMP & SIG AS THEY ARE FITTED C VALUES OUT OF LEAST SQUARES ANALYSIS C C SPECIAL NOTE THIS ANALYSIS IS BASED ON THE FOLLOWING CONCEPT C C IF Y(EI) = GAUSSIAN (EI) THEN C LN (Y) = LN (AMP) - (E-EBAR)**2/(2*SIG**2) C THIS IS OF THE FORM C Y = B + X * M C IF YOU IDENTIFY B = LN(AMP) AND M= -1/(2*SIG**2) C AND YOU ASSUME THAT THE INDEPENDENT VARIABLE X = (E-EBAR)**2 C C IT WILL ONLY WORK WELL IF THE INITIAL ESTIMATE OF EBAR IS C ESSENTIALLY THE TRUE GAUSSIAN CENTROID OF THE PEAK AND ALSO C IF THE PEAK IS A SYMMETRIC GAUSSIAN, IF THERE IS C SIGNIFICANT LOW ENERGY TAILING OF THE GAUSSIAN DUE TO C INCOMPLETE CHARGE COLLECTION AND/OR SUM PEAKS THEN THE FIT C WILL LIKELY BE POOR C C DO 1000 I = IBL,IBU X=DBLE((BX(I)-EBAR)**2) Y=DLOG(DBLE(BY(I))) CALL LSTSQR(2,X,Y,,,,SUM) 1000 CONTINUE CALL LSTSQR(3,,,NUMPTS,sM,B,SUM) AMP=EXP(B) SIG=1./SQRT(-2.*sM) RETURN END C----------------------------------------------------------------------- C C SUBROUTINE LSTSQR(ICOM,X,Y,NUMPTS,SLOPE,YINTER,SUM) C SOURCE: NESTOR ZALUZEC C$ FILENAME: ELLTSQ.FOR C$ PROJECT: PV9100/60 C$ SUB-PROJECT: EELS C$ VERSION: 1.0 C$ ENVIRONMENT: 2.2 C$ CREATED: 15-OCT-81 BY W. Bresler C$ EDIT LOG: C$ NO. DATE BY REASON C$ 1 25-AUG-82 WB Use double precision for better accuracy C 2 15-FEB-83 NJZ FOR USE ON GENERAL RT-11 ROUTINE WITH NO C CLRA FNCT C$ DOCUMENTATION: C C ********************************************************************** C * C * This routine will perform least squares fitting C * to a linear function. C * C * Parameters: C * Inputs: C * ICOM - the command parameter: C * 1- zeros the summations C * 2- adds a point to the summations C * 3- computes the slope and y-intercept C * 4- computes the correlation coeff. C * NUMPTS - the number of points in the fit. C * needed only when ICOM = 3,4 C * X - the value of x in an (x,y) pair C * Y - the value of y in an (x,y) pair C * SLOPE - the returned value input to compute C * the corr. coeff. C * SUM - the array used to hold the summations C * C * Outputs: C * SLOPE - the slope of the fit; returned when C * ICOM = 3. C * YINTER - the y-intercept of the fit returned C * when ICOM = 3; C * the corr. coeff of the fit returned C * when ICOM = 4 C * C * Errors: C * NUMPTS = 1 C * Result - SLOPE = 0 ; YINTER = y C * Action - none necessary C * C ********************************************************************** C C DOUBLE PRECISION SUM(5),X,Y,DENOM,SIGMAX,SIGMAY C C C Zero Sigma's Fit Corr coeff GO TO (1100, 1110, 1120, 1130)ICOM C ------------------------------------------------------------ C ZERO ALL SUMMATIONS C ------------------------------------------------------------ 1100 DO 1101 I=1,5 1101 SUM(I)=0. GO TO 1150 !Return point C ------------------------------------------------------------ C ADD A POINT TO THE FIT C ------------------------------------------------------------ 1110 SUM(1) = SUM(1) + X !SIGMA(X'S) SUM(2) = SUM(2) + Y !SIGMA(Y'S) SUM(3) = SUM(3) + X * X !SIGMA(X*X) SUM(4) = SUM(4) + X * Y !SIGMA(X*Y) SUM(5) = SUM(5) + Y * Y !SIGMA(Y*Y) GO TO 1150 C ------------------------------------------------------------ C COMPUTE THE FIT PARAMETERS C ------------------------------------------------------------ 1120 IF ( NUMPTS .GT. 1 ) GOTO 1125 SLOPE = 0 ! Single point fit YINTER = SUM(2) GO TO 1150 1125 DENOM = NUMPTS * SUM(3) - SUM(1) * SUM(1) ! Multiple point fit SLOPE = (SUM(4) * NUMPTS - SUM(1) * SUM(2))/DENOM YINTER = (SUM(3) * SUM(2) - SUM(1) * SUM(4))/DENOM GO TO 1150 C ------------------------------------------------------------ C COMPUTE THE CORRELATION COEFFICIENT C ------------------------------------------------------------ 1130 SIGMAX = DSQRT( SUM(3) / NUMPTS - ( SUM(1) / NUMPTS ) ** 2 ) SIGMAY = DSQRT( SUM(5) / NUMPTS - ( SUM(2) / NUMPTS ) ** 2 ) YINTER = SLOPE * SIGMAX / SIGMAY !RETURNS VALUE IN YINTER GO TO 1150 C ------------------------------------------------------------ C NORMAL RETURN PATH C ------------------------------------------------------------ 1150 RETURN END C ___________________________________________________________ SUBROUTINE SHIFT(A,MM,ICENT) C __________________________________________________________ C C SHIFTS AN ARRAY SO THAT ICENT IS THE FIRST CHANNEL. C ALL CHANNEL LESS THAN ICENT ARE SHIFTED TO THE LAST CHANNELS C THUS THERE IS WRAP AROUND AND NO DATA LOST. C MM - SIZE OF ARRAY C A - ARRAY TO BE SHIFTED C ICENT - CHANNEL IN ARRAY TO BE FIRST CHANNEL DIMENSION A(4096),STORE(4096) DO 100 I=1,ICENT-1 100 STORE(I+MM-ICENT+1)=A(I) DO 200 I=ICENT,MM 200 STORE(I-ICENT+1)=A(I) 250 FORMAT(I5,2F10.4) DO 300 I=1,MM 300 A(I)=STORE(I) RETURN END C C C ********************************************************* C SUBROUTINE ZEROFIT(MMIN,NN,NZ,AMP,SIGFIT,IFWTM,IBL) C C THIS SUBROUTINE FITS A GAUSSIAN FUNCTION TO THE ZERO LOSS C PEAK BY THE LEAST SQUARES METHOD. C C THIS ROUTINE CALLS ON THE SUBROUTINE ZLPKFIT FOR A C LEAST SQARES FIT OF THE ZEROLOSS PEAK TO A GAUSSIAN FORM. C THE VALUE OF SIGMA FROM THE LEAST SQUARES FIT IS SIGFIT. THE USER C IS INFORMED OF THE CALCULATED VALUE OF THE SIGMA OF THE ZERO LOSS C PEAK (SIGFIT) AND IS REQUESTED TO INPUT A SIGMA TO BE USED AS THE C SIGMA OF THE UNIT AREA GAUSSIAN USED IN THE DECONVOLUTION C COMMON /D/ D DIMENSION BX(4096),BY(4096) DIMENSION D(8192) C C CALLS ON SUBROUTINE ZLPKFIT TO FIT ZEROLOSS PEAK TO A GAUSSIAN C SET UP FOR SUBROUTINE ZLPKFIT: C C CREATE X AXIS DO 600 I=1,NN 600 BX(I)=FLOAT(I) C C DETERMINE UPPER AND LOWER BOUNDS OF RANGE OF ZERO LOSS C ARRAY TO BE FIT BY ZLPKFIT ROUTINE. C THE FULL WIDTH TENTH MAXIMUM OF EITHER SIDE OF THE C ZERO LOSS PEAK SETS THE UPPER AND LOWER BOUNDS FOR THE FIT. FWTM=D(2*NZ-1)/10. DO 601 I=1,NZ IF(D(2*I-1).GE.FWTM) GOTO 602 601 CONTINUE 602 IBL=I DO 603 I=NZ,NN IF(D(2*I-1).LE.FWTM) GOTO 604 IFWTM=I 603 CONTINUE 604 IBU=I-1 IF(IFWTM.GT.MMIN) THEN IFWTM=MMIN IBU=IFWTM-1 ENDIF WRITE(6,610) 610 FORMAT(' HIT TO CONTINUE ') READ(6,615)X 615 FORMAT(1A80) WRITE(6,606) 606 FORMAT('0ZERO LOSS PEAK STATISTICS:',/) WRITE(6,607)IBL,IBU 607 FORMAT('0LOWER AND UPPER BOUND CHANNELS USED FOR ',/,'0FITTING @ THE ZERO LOSS PEAK ',2I5) EBAR=FLOAT(NZ) C C CREATE Y AXIS C DO 605 I=IBL,IBU 605 BY(I)=D((2*I)-1) C C C CALL ZLPKFIT(BX,BY,IBL,IBU,AMP,EBAR,SIGFIT) C C RETURN END C----------------------------------------------------------------------- SUBROUTINE RPDL(ND,EPC,E0,BETA,IOUTNAME) C SOURCE: NESTOR ZALUZEC C----------------------------------------------------------------------- C SUBROUTINE TO READ AN ASCII DATA FILE IN EMMPDL DATA FORMAT C INTO AN ARRAY XS(4110) C C FIRST TWO LINES OF OF FILE CONTAIN TEXT ARE FOR IDENTIFICATION C AND CAN BE READ USING A 80A1 FORMAT C THE NEXT LINE CONTAINS THE EXPERIMENTAL DATA STORED AS 10 ASCII C NUMBERS (FLOATING POINT FORMAT) EACH DATA POINT SEPERATED BY A C COMMA AS A DELIMINATOR, C THE FIRST NUMBER OF THE FIRST LINE IS THE TOTAL NUMBER OF DATA C POINTS IN THE SPECTRAL FILE, AND IS LIMITED TO A MAXIMUM VALUE C OF 4095, IF IT'S VALUE IS "ND" THEN THE MAXIMUM NUMBER OF C LINES OF DATA IN THIS FILE IS GIVEN BY NLINES="ND/10" +1 C FOLLOWING NLINES OF DATA THERE IS A FINAL LINE OF 10 VALUES C WHICH ARE THE CALIBRATION CONSTANTS OF THE SPECTRUM THEY C ARE DEFINED AS FOLLOWS: C C XS(1) = ND = NUMBER OF POINTS (CHANNELS) IN THE SPECTRUM C XS(2)-XS(ND) =DATA C C XS(ND+1) TO XS(ND+10) =CALIBRATION CONSTANTS OF SPECTRAL DATA C ND+1= OFFSET ENERGY OF FIRST CHANNEL C ND+2= EV/CHANNEL C ND+3= ACCELERATING VOLTAGE (KV) C ND+4= INCIDENT BEAM DIVERGENCE (MR) FOR EELS DATA C = X-AXIS (ROD AXIS) TILT (DEGREES) FOR XEDS DATA C ND+5= SCATTERING ANGLE (MR) (BETA) C = Y-AXIS (CUP AXIS) TILT (DEGREES) FOR XEDS DATA C ND+6= LIVE TIME (DWELL*NUMBER OF CHANNELS) FOR EELS OR C SIMPLE LIVE TIME FOR XEDS DATA C ND+7= DEAD TIME C ND+8= INCIDENT BEAM CURRENT IN nA C ND+9= PROBE DIAMETER IN nm C ND+10= SPECIMEN THICKNESS ALONG BEAM DIRECTION IN nm. C C C C----------------------------------------------------------------------------- C COMMON /RPDL/XS DIMENSION INAME(7) DIMENSION XS(4110),IXL(80) BYTE DAY(9),TIM(8) C C READ IN DATA FILE NAME C 204 WRITE(6,1) 1 FORMAT(' ENTER DATA FILE NAME (FILENAME.DAT):',/) READ(6,205) (INAME(J),J=1,7) 205 FORMAT(7A2) C OPEN (UNIT=11,NAME=INAME,TYPE='OLD',ACCESS='SEQUENTIAL') C C C WRITE IT OUT C C CALL DATE(DAY) CALL TIME(TIM) WRITE(6,206) DAY,TIM,(INAME(J),J=1,7) 206 FORMAT(1X,9A1,' at ',8A1,/,' Input Data File = ',7A2) C C READ IN TWO LINES OF TEXT C C C READ(11,2) IXL 2 FORMAT(80A1) WRITE(6,3) IXL 3 FORMAT(1X,80A1) READ(11,2) IXL WRITE(6,3) IXL C C READ FIRST DATA LINE AND DETERMINE THE NUMBER OF CHANNELS C READ (11,1000) (XS(I), I=1,10) 1000 FORMAT (10F12.0) ND=XS(1) WRITE (6,1001) ND 1001 FORMAT (/,' THIS SPECTRAL FILE CONTAINS ',I4,' CHANNELS @ OF DATA') C C CALCULATE THE REMAINING NUMBER OF LINES (YOU HAVE ALREADY READ ONE!) C IF((XS(1)/10. - ND/10 ).GT.0) NLINES= ND/10 +1 IF((XS(1)/10. - ND/10 ).EQ.0) NLINES= ND/10 C C NOW LOOP THROUGH THE DATA READING EACH LINE C DO 4 I=2,NLINES IF=10*I-9 IL=10*I IF (IL.GE.ND) IL=ND C READ(11,1000,END=4,ERR=4) (XS(J), J=IF,IL) 4 CONTINUE C C READ IN CALIBRATION CONSTANTS C IF THEY EXIST C READ(11,1000,ERR=52,END =52) (XS(J), J =ND+1,ND+10) 52 CONTINUE IF((XS(ND+1).EQ.0).AND.(XS(ND+2).EQ.0)) WRITE(6,53) 53 FORMAT(' UNCALIBRATED SPECTRUM:') WRITE(6,601) 601 FORMAT (/,' EXPERIMENTAL CALIBRATION CONSTANTS: ') WRITE(6,7) (XS(J), J=ND+1,ND+5) 7 FORMAT(//, C ' OFFSET = ',1F9.3,' eV',/, C ' EV/CHANNEL = ',1F9.3,/, C ' ACCELERATING VOLTAGE = ',1F9.1,' kV',/, C ' INCID. BEAM DIVERGENCE = ',1F9.3,' mR',/, C ' SCATTERING ANGLE = ',1F9.3,' mR',/, C ' LIVE TIME = ',1F9.3,' SEC',/, C ' DEAD TIME = ',1F9.3,' SEC',/, C ' BEAM CURRENT = ',1F9.3,' nA',/, C ' PROBE DIAMETER = ',1F9.3,' nm',/, C ' SPECIMEN THICKNESS = ',1F9.3,' nm',/) 6 WRITE(6,600) 600 FORMAT (' HIT TO CONTINUE') READ(6,2) X CALL CLOSE(11) EPC=XS(ND+2) E0=XS(ND+3) BETA=XS(ND+5) C RETURN END C C ************************************************************* C SUBROUTINE SPECEXT(MFIN,A,R) C THIS SUBROUTINE DETERMINES EXTRAPOLATION VALUES A AND R C FOR EXTRAPOLATING SPECTRA OF THE FORM A*E**R COMMON/D/ D DIMENSION D(8192) C REAL*8 X(100),Y(100) DOUBLE PRECISION SUM(5),X,Y C C FIT LAST TEN DATA POINTS NUMPTS=10 C C CREATE X AXIS C PLACE NATURAL LOG OF DATA IN Y ARRAY C C CALL LEAST SQUARES ROUTINE TO FIT LINEARIZED DATA C OF THE FORM LN(D)=LN(A)+R*LN(E) C VALUE RETURNED FOR SLOPE WILL BE THE VALUE FOR R C VALUE RETURNED FOR INTERCEPT WILL BE VALUE FOR LN(A) CALL LSTSQR(1,,,,,,SUM) DO 100 I=1,NUMPTS X=DLOG(DFLOAT(I+(MFIN+1)/2-NUMPTS)) Y=DLOG(DBLE(D(2*I+MFIN-2*NUMPTS))) CALL LSTSQR(2,X,Y,,,,SUM) 100 CONTINUE CALL LSTSQR(3,,,NUMPTS,SLOPE,YINTER,SUM) R=SLOPE WRITE(6,198) 198 FORMAT(//,' THE END OF THE SPECTRUM WILL BE FITTED TO A FUNCTION @ OF THE FORM A*EXP(-R),') WRITE(6,199) 199 FORMAT (' WHERE VALUES OF A AND R ARE FITTED ACCORDING @ TO THE DATA.',/,' R SHOULD BE -5 < R < -2.',//) WRITE(6,200)R 200 FORMAT(' RESULTS FROM FITTING THE SPECTRUM GAVE A VALUE OF R= @',F10.3) C 205 if(r.gt.-2. .or. r.lt.-5.) then WRITE(6,210)R 210 FORMAT('1THE EVALUATED VALUE FOR R,',F7.1,' IS NOT REASONABLE') WRITE(6,211) 211 FORMAT(' PLEASE ENTER VALUE FOR R BETWEEN -2. AND -5.') READ(6,*)R A=D(MFIN)/(((MFIN+1)/2)**R) WRITE(6,200)R if(r.gt.-2. .or. r.lt.-5.) GO TO 205 ELSE A=EXP(YINTER) GO TO 300 ENDIF 300 RETURN END C ************************************************************* C C ___________SAMPLE DATA FILE "AL.DAT" STARTS AFTER THIS LINE____________ Sample: Aluminum Spectrum = Example from the text "EELS in the E.M." by R.Egerton 79, 2, 4, 7, 25, 166, 489, 702, 539, 261, 119, 61, 38, 31, 31, 31, 32, 33, 37, 50, 103, 236, 420, 450, 329, 196, 129, 94, 77, 70, 67, 64, 60, 59, 59, 67, 92, 144, 179, 178, 137, 103, 80, 68, 61, 57, 55, 51, 49, 47, 45, 46, 52, 61, 68, 65, 58, 49, 42, 38, 36, 34, 32, 31, 29, 28, 27, 27, 26, 28, 28, 27, 24, 22, 21, 20, 19, 18, 18, 17, 0,1., 100, 5.3, 20, 0, 0, 1000,