write_mupho.pro 12.9 KB
Newer Older
1
FUNCTION write_mupho, w_in, file=file, phimin=phimin, phimax=phimax, abs=abs, $
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
  aemp=aemp, emax = emax, sam=sam, npho=npho, itm=itm, $
  ivit=ivit, idw=idw, iemp=iemp, ires=ires, ipr=ipr, iloss=iloss, $
  amass=amass, sig=sig, conc=conc, $
  unt=unt, fun=fun
  ;+
  ;
  ; write_mupho(w_in, file='input_mupho.txt')
  ;
  ; w_in: Spectrum in time-of-flight
  ;
  ; write_mupho sum the sin(theta)-weighted intensities of the instrument-free
  ; time-of-flight data and create the input_mupho.txt input file for the routine
  ; MUPHOCOR.
  ;
  ; prox example:
  ;
  ; RDSET,inst="IN5"
  ; RDSET,base="Current Path"
  ; w_clear, /all
  ;
  ;  ;  read the time-of-flight reduced data
  ;  ; BGA at room temperature
  ;  w1 = rdrun('152965_152990.hdf')
  ;
  ;  ; Create the MUPHOCOR input with various options
  ;  w2 = write_mupho(w1, file='input_mupho.txt', emax=40.0, abs=0.0, phimin=31.0, phimax=134.8, iemp=0, amass=44.8, sig=3.9, conc=1.0, unt=0.015)
  ;
  ;  ; Compute the GDOS with the MUPHOCOR routine
  ;  w10 = muphocor('input_test.txt') & see, w10, /below
  ;
  ;
  ;
  ; input_mupho.txt example and parameters:
  ;
  ;!     E0:INC. ENERGY(MEV)             FP:FLIGHTPATH(CM)
  ;!     XNEL:CHANNEL OF ELASTIC LINE    CW:CHANNEL WIDTH(MYS)
  ;!     FIMI:MIN.SCATT.ANGLE            FIMA:MAX.SCATT.ANGLE
  ;!     UNT:CONST.BACKGROUND            ABK: Absorption coefficient
  ;!     AEMP: Coeff. for calculating the detector efficiency
  ;
  ; e0,     fp,         cw,      xnel,     fimi,    fima,       abk,      aemp
  ; 4.77    248.30      3.91    794.08     10.33    115.05      0.00      3.30
  ;
  ;
  ;!       Correction OF constant backgroud FOR counter efficiency
  ;  bac=unt/(1.0-EXP(-aemp/SQRT(e0-hot(n))))
  ;
  ;
  ;
  ;!     NSPEC: NUMBER OF DIFFERENT ATOMIC SPECIES
  ;!     IF IQUO is unequal 0 an external model FOR the ratio OF F/G has
  ;!     to be provided. The model is given through the array PARQ
  ; lm,     nspec,  iquo,lpar
  ; 500     1       0    1
  ;
  ;!     HOX:   UPPER LIMIT OF DENSITY OF STATES
  ;!     TEMPO: TEMPERATURE IN KELVIN
  ;!     DW:    MEAN DEBY WALLER COEFFICIENT
  ;
  ;hox,    temp0,      dw
  ;40.00    298.52      0.03
  ;
  ;
  ;!     AMASI: ATOMIC MASS
  ;!     SIGI : SIGMA
  ;!     CONCI: CONCENTRATION
  ;!     ALFI:  SCATTERING POWER
  ;DO  i=1,nspec <-- 1 in our CASE so far
  ;READ 40,amasi(i),sigi(i),conci(i)
  ;END DO
  ;amasi(i),  sigi(i),  conci(i)
  ;44.80      3.90      1.00
  ;
  ;
  ;!     NPHO:     NUMBER OF MULTI PHONON TERMS
  ;!     ITM:      TOTAL NUMBER OF ITERATIONS
  ;!     IVIT=0(1):ITERATION BY DIFFERENCE(QUOTIENT) METHOD
  ;!     IDW=0:    DW COEFF. KEPT CONST.=INPUT VALUE,DW=1:DW COEFF.ITERATED
  ;!     IEMP=1(0):CORRECTIONS FOR COUNTER EFFICIENCY WILL BE(NOT BE) DONE
  ;!     IRES=1(0):DATA WILL BE(NOT BE) CORRECTED FOR SPECTR. RESOLUTION
  ;!     IPR=1(0): Convolutions integrals o the multiphonon term are printed
  ;!     ILOSS=1:  Analysis OF the imcomplete energy loss spectrum will be done
  ;
  ;npho,itm, ivit, idw, iemp, ires, ipr,iloss
  ;5    10   10    1    1     0     0   0
  ;
  ;
  ;
  ;!     NU:  First channel OF spektrum NO: Last channel OF spektrum
  ;!     NUU: First channel OF TOF distribution used FOR calculation
  ;!     NOO: Last channel  OF TOF distribution used FOR calculation
  ;!     IGLU: Number OF smoothing processes in CASE OF a time-dependant background
  ;
  ;isour,nuu, noo, nu,  no,    iglu
  ;0     200  750  200  750    0
  ;
  ;
  ;!     UNT: CONSTANT BACKGROUND
  ;!     FUN: MULTIPLICATION FACTOR FOR TIME DEPENDENT BACKGROUND
  ;
  ;unt,      fun
  ;1.0000    0.0000
  ;
  ;
  ;;     Z00: TIME OF FLIGHT DISTRIBUTION
  ;;     UN0: TIME DEPENDENT BACKGROUND (if fun = 1)
  ;

  ; :Author: ollivier (2016-12-23)
  ;-

  ;Corrections:
  ;
  ;2017-07-04: doesn't work for water: no elastic line large enough

117
118


119
120

  COMMON c_lamp_access, inst
121
122

  s = SIZE(w_in)
123
124
  IF s[0] NE 2 THEN GOTO, ERROR

125
126
127
128
129
  TAKE_DATP, datp
  par  = datp.p
  x_in = datp.x
  y_in = datp.y
  title= datp.w_tit
130
131
132


  IF (NOT KEYWORD_SET(phimin)) THEN phimin = min(y_in)
133
134
  IF (NOT KEYWORD_SET(phimax)) THEN phimax = max(y_in)
  idx = where(y_in GT phimin AND y_in LT phimax)
135

136
137
  phimin = y_in[idx[0]]
  phimax = y_in[idx[n_elements(idx)-1]]
138

139
140
141
  ; ----------------------------------------------------------------------------
  ; From the Reichardt report:
  ; --------------------------
142
143
  ; "The analysis performed by MUPHCOR starts from the sin(theta) weighted
  ; angle integrated time-of-flight distribution."
144
145
  ; WARNING: y_in is 2theta => theta = y_in/2
  ; ----------------------------------------------------------------------------
146
147
  ;  w_buf = w_in[*,idx] * ((fltarr(s[1])+1.0)#y_in[idx]/2.0/180.0*!PI)
  w_buf = w_in[*,idx] * sin((fltarr(s[1])+1.0)#y_in[idx]/2.0/180.0*!PI)
148
  w_buf = total(w_buf,2)
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175

  ; ----------------------------------------------------------------------------
  ;  Not to do?: the position of ELP gives the information to the routine
  ;  From other tests input this is not to be done.
  ;  IF STRUPCASE(inst) NE 'IN4' THEN  w_buf = REVERSE(w_buf)
  ;
  ;  2017-07-04
  ;  No way to check from within how long distance from the ELP to take data
  ;  to improve urgently
  ; ----------------------------------------------------------------------------
  idx = (par[9] - 30)>0.0 + findgen(60)
     
  gs = GAUSSFIT(x_in[idx], w_buf[idx], coeff, NTERMS=3)

;Usual case, high T: use the AS side
   idx = where(x_in LE coeff[1]-coeff[2]/2)

;
;  print, '****************************************************************'
;  print, 'WARNING: Stokes case, see in  
;  print, '/home/cs/lambda/macros/IN5/JOR_LIBRARY/write_mupho.pro'
;  print, 'line `160 and folowing to reverse to normal situation (Anti-Stokes)'
;  print, '****************************************************************'
;
;Unusual case, low T: use the Stokes side 
;-----------------------------------------
;  idx = where(x_in GE coeff[1] + 2.5*coeff[2])
176
  
177
;  print, coeff(1) + 2*coeff(2)
178
  
179
180
  imin = idx[0]
  imax = idx[n_elements(idx)-1]
181
  
182
  print, 'Min & max index for calculation: ', imin, imax
183

184
185
186
187
188
189
190
;Do it twice otherwise it does not do anything ???...
  w_buf = w_buf[idx]
;  w_buf = w_buf[idx]
 
; print, 'Min & max index for calculation: ', imin, imax
 
 
191
192
193
194
195
  ; ---------------------------------------------------------------------------
  ;     E0:INC. ENERGY(MEV)             FP:FLIGHTPATH(CM)
  ;     XNEL:CHANNEL OF ELASTIC LINE    CW:CHANNEL WIDTH(MYS)
  ;     FIMI:MIN.SCATT.ANGLE            FIMA:MAX.SCATT.ANGLE
  ;     UNT:CONST.BACKGROUND            ABK: Absorption coeff. (0.0 ?)
196
  ;     AEMP: Coeff. for calculating detector efficiency:
197
198
199
200
  ;           Correction of constant backgroud for counter efficiency
  ;     IF(iemp == 0) bac=unt/(1.0-EXP(-aemp/SQRT(e0-hot(n))))
  ;     AEMP = 3.3 in De Boissieu data and also in others fo IN6 data
  ; ---------------------------------------------------------------------------
201
  E0=0. & FP=0. & CW=0. & XNEL=0. & FIMI=0. & FIMA=0. & ABK=0.
202
203
204
  E0 = 81.8045/par[21]^2   ; incident energy
  FP = par[27]*100         ; unit [cm]
  CW   = par[18]           ; channel width [mu-s]
205
  XNEL = par[9]            ; ELP (as computed after "normalise")
206
207
208
209
  FIMI = phimin            ; min. scatering angle
  FIMA = phimax            ; max. scatering angle
  IF (NOT KEYWORD_SET(abs)) THEN abs = 0.
  ABK  = abs
210
211

  ; Coef for detector efficiency: value from De Boissieu et al. taken (IN6 data)
212
213
  IF (NOT KEYWORD_SET(aemp)) THEN AEMP = 3.3

214

215
216
217
218
219
  ; ---------------------------------------------------------------------------
  ;     NSPEC: NUMBER OF DIFFERENT ATOMIC SPECIES
  ;     LM: Number of points
  ;     IQUO: Indicator to select model F(w)/G(w)  = 0
  ;     LPAR: Nb parameters for calculation of  F(w)/G(w) = 1
220
  ;
221
222
223
224
  ;     LM must be limited to 512 (LM2 and LM5 to 1024 in MUPHOCOR)
  ; ---------------------------------------------------------------------------
  LM = n_elements(w_buf)<512

225
226
227
  print, '-------------------------------------------------------'
  print, 'Number of points LM = ', strtrim(string(LM),2)
  print, '-------------------------------------------------------'
228
229
230
231
232
233
234
235
236
237
238
239
240

  NSPEC=1L & IQUO=0L & LPAR= 1L

  ; ---------------------------------------------------------------------------
  ;     HOX:   UPPER LIMIT OF DENSITY OF STATES
  ;     TEMPO: TEMPERATURE IN KELVIN
  ;     DW:    MEAN DEBY WALLER COEFFICIENT
  ; ---------------------------------------------------------------------------
  HOX=0. & TEMP0=0. & DW=0.

  IF (NOT KEYWORD_SET(emax)) THEN HOX = 100. ELSE HOX = emax
  IF (NOT KEYWORD_SET(DW)) THEN DW = 0.
  IF par[11] GT 0. THEN TEMP0 = par[11] ELSE TEMP0 = 295.
241

242
243
244
245
246
247
248
249
250
  ; ---------------------------------------------------------------------------
  ;     AMASI: ATOMIC MASS
  ;     SIGI : SIGMA
  ;     CONC:  CONCENTRATION
  ;     ALFI:  SCATTERING POWER
  ; ---------------------------------------------------------------------------
  SIGI =fltarr(2) & CONCI=SIGI & AMASI=SIGI
  IF (NOT KEYWORD_SET(emax)) THEN emax=max(x_in)
  sam = fltarr(3)
251
  IF (NOT KEYWORD_SET(amass)) THEN amass = 123.0
252
253
254
  IF (NOT KEYWORD_SET(sig)) THEN sig = 3.9
  IF (NOT KEYWORD_SET(conc)) THEN conc = 1.0
  IF (NOT KEYWORD_SET(amass)) OR  (NOT KEYWORD_SET(sig)) OR (NOT KEYWORD_SET(conc)) THEN $
255
256
    print, 'WARNING: Sample parameters: Sigma, concentration or Atomic mass not defined!'

257
258
259
260
261
262
263
  sam[0] = amass & sam[1] = sig & sam[2] = conc

  FOR I=0,NSPEC-1 DO BEGIN
    AMASI[I] = sam[0]
    SIGI[I]  = sam[1]
    CONCI[I] = sam[2]
  ENDFOR
264
265


266
267
268
269
270
271
272
273
  ; ---------------------------------------------------------------------------
  ;     NPHO:     NUMBER OF MULTI PHONON TERMS,ITM:NUMBER OF ITERATIONS
  ;     ITM:      TOTAL NUMBER OF ITERATIONS
  ;     IVIT=0(n):ITERATION BY DIFFERENCE(QUOTIENT) METHOD UP TO iteration n
  ;     IDW=0:    DW COEFF. KEPT CONST.=INPUT VALUE,DW=1:DW COEFF.ITERATED
  ;     IEMP=1(0):CORRECTIONS FOR COUNTER EFFICIENCY WILL BE(NOT BE) DONE
  ;     IRES=1(0):DATA WILL BE(NOT BE) CORRECTED FOR SPECTR. RESOLUTION
  ; ---------------------------------------------------------------------------
274
  IF (NOT KEYWORD_SET(nphon)) THEN  NPHO = 5L
275
276
277
278
279
280
281
  IF (NOT KEYWORD_SET(itm))   THEN  ITM  = 10L
  IF (NOT KEYWORD_SET(ivit))  THEN  IVIT = 10L
  IF (NOT KEYWORD_SET(idw))   THEN  IDW  = 1L
  IF (NOT KEYWORD_SET(iemp))  THEN  IEMP = 0L
  IF (NOT KEYWORD_SET(ires))  THEN  IRES = 0L
  IF (NOT KEYWORD_SET(ipr))   THEN  IPR  = 0L
  IF (NOT KEYWORD_SET(iloss)) THEN  ILOSS= 0L
282

283
284
285
286
287
288
  ; ---------------------------------------------------------------------------
  ;     NU:  First channel of spektrum NO: Last channel of spektrum
  ;     NUU: First channel of TOF distribution used for calculation
  ;     NOO: Last channel  of TOF distribution used for calculation
  ;
  ; ---------------------------------------------------------------------------
289
290
291
292
293
294
295
296
297
298
299
300
301
  ISOUR=0L
  ;  MUPHOCOR doesn't work well with neutron energy loss (ENERGY LOSS mode)???
;    IF (NOT KEYWORD_SET(nuu)) THEN  $
;      IF STRUPCASE(inst) EQ 'IN4' THEN NUU=imax[0]+dN+1 ELSE NUU = 1L
;    IF (NOT KEYWORD_SET(noo)) THEN  $
;      IF STRUPCASE(inst) EQ 'IN4' THEN NOO= n_elements(w_buf)-1+NUU ELSE NOO=n_elements(w_buf)

;  NUU = imin
;  NOO = imax

  ; ---------------------------------------------------------------------------
  ; this is what we have usally (in case AS)
  ; ---------------------------------------------------------------------------
302
  NUU = 1L
303
  NOO = n_elements(w_buf)
304
305
  NU=NUU
  NO=NOO
306
  ; ---------------------------------------------------------------------------
307
  IGLU=0L
308

309
310
311
312
313
314
  ; ---------------------------------------------------------------------------
  ;     UNT: CONSTANT BACKGROUND
  ;     FUN: MULTIPLICATION FACTOR FOR TIME DEPENDENT BACKGROUND
  ; ---------------------------------------------------------------------------
  IF (NOT KEYWORD_SET(unt)) THEN  UNT=0.0
  IF (NOT KEYWORD_SET(fun)) THEN  FUN=0.0
315

316
317
318
319
  ; ---------------------------------------------------------------------------
  ;     Z00: TIME OF FLIGHT DISTRIBUTION (spectrum 1)
  ;     UN0: TIME DEPENDENT BACKGROUND   (spectrum 2?)
  ; ---------------------------------------------------------------------------
320
321
    Z00=fltarr(1024) & UN0=Z00

322
323
324
325
326
327
328
329
330
331
332
333
334
335
  ; ---------------------------------------------------------------------------
  ; open the output file and write parameters and intensities
  ; ---------------------------------------------------------------------------
  OPENW,  lun, file, /GET_LUN
  PRINTF, lun, title,format = "(A)"
  PRINTF, lun, E0, FP, CW, XNEL, FIMI, FIMA, ABK, AEMP, format="(8F10.2)"
  PRINTF, lun, LM, NSPEC, IQUO, LPAR, format="(4I5)"
  PRINTF, lun, HOX, TEMP0, DW, format="(3F10.2)"
  FOR I=0,NSPEC-1 DO BEGIN
    PRINTF, lun, AMASI[I], SIGI[I], CONCI[I], format="(3F10.2)"
  ENDFOR
  PRINTF, lun, NPHO,ITM,IVIT,IDW,IEMP,IRES,IPR,ILOSS, format="(8I5)"
  PRINTF, lun, ISOUR,NUU,NOO,NU,NO,IGLU, format="(6I5)"
  PRINTF, lun, UNT, FUN, format="(2F10.4)"
336

337
338
339
340
341
342
343
344
345
  IF (ISOUR NE 1) THEN BEGIN
    PRINTF,lun, w_buf, format="(8F10.4)"
    IF (FUN NE 0.) THEN BEGIN
      PRINTF,lun,tmp,format="(2F10.4)"
    ENDIF
  ENDIF
  PRINTF, lun, 0., format="(1F10.4)"


346
  ERROR:
347
  FREE_LUN, lun
348
349
350
  
;  In all cases, return
  RETURN, w_buf
351
352

END