Документ взят из кэша поисковой машины. Адрес оригинального документа : http://hea-www.harvard.edu/PINTofALE/pro/loopem.pro Дата изменения: Tue Jun 28 01:52:23 2005 Дата индексирования: Tue Oct 2 00:13:27 2012 Кодировка: Поисковые слова: annular solar eclipse

function loopem,logTmax,EM,logT=logT,sloop=sloop,verbose=verbose, _extra=e
;+
;function loopem
; convert a set of integrated emission measures to a differential
; emission measure distribution assuming a hydrostatic-loop type
; distribution and return computed DEM
;
; DEM(T) is generally n^2*dX/dlnT, and
; EM(T) == \int_0^T dlnT' DEM(T') = \int_0^T dT' DEM(T')/T'
; In particular, DEM is parameterized as DEM(T)=a*T^b
; So for a given EM(T), we obtain
; EM(T)=(a/b)*(Tmax^(b)-Tmin^(b)), Tmin=0, b>0
; with Tmin=0 for b>0
;
;syntax
; dem=loopem(logTmax,EM,logT=logT,sloop=sloop,verbose=verbose)
;
;parameters
; logTmax [INPUT; required] loop top temperature(s) in log10(deg K)
; EM [INPUT] if given, normalize the output DEM to give a
; total Emission Measure equal to this.
; * size must match that of logTmax -- if fewer, fills out
; with 1st element; if more, ignored
;
;keywords
; logT [INPUT] temperature grid over which output DEM is returned
; * default is 4..8 in steps of 0.05
; sloop [INPUT] slope of DEM in log(DEM)-log(T) space
; * default is 1.5
; verbose [INPUT] controls chatter
; _extra [JUNK] here only to prevent crashing the program
;
;history
; vinay kashyap (Nov01)
; changed NORM calc to allow for large values of SLOOP (VK; Feb05)
;-

; usage
ok='ok' & np=n_params() & ntmx=n_elements(logTmax)
if np eq 0 then ok='Insufficient parameters' else $
if ntmx eq 0 then ok='logTmax is undefined'
if ok ne 'ok' then begin
print,'Usage: dem=loopem(logTmax,EM,logT=logT,sloop=sloop,verbose=verbose)'
print,' return DEM(T) [cm^-5/logK] of a static RTV loop'
if np ne 0 then message,ok,/info
return,-1L
endif

; check inputs
vv=0 & if keyword_set(verbose) then vv=long(verbose[0]) > 1
;
logTm=logTmax
;
EMtot=dblarr(ntmx)+1.0 & nem=n_elements(EM)
if nem gt 0 then begin
EMtot=dblarr(ntmx)+double(EM[0])
if nem le ntmx then EMtot[0L:nem-1L]=EM else EMtot=EM[0L:ntmx-1L]
endif
;
tlog=findgen(81)*0.05+4. & nT=n_elements(logT)
if nT ne 0 then tlog=logT & nT=n_elements(tlog)
DEM=dblarr(nT)
;
ss=1.5 & if keyword_set(sloop) then ss=sloop[0]
ok='ok'
if ss le 0 then begin
message,'slope cannot be -ve',/info
if vv gt 4 then message,'Assuming delta-fn T distributions',/info
for i=0L,ntmx-1L do begin
jnk=max(abs(logTm[i]-tlog),imx)
dlogT=(logT[imx]-logT[[imx-1L]]) > (logT[[imx+1L]]-logT[imx])
if dlogT eq 0 then dlogT=1.
DEM[imx]=DEM[imx]+EMtot[i]/dlogT
endfor
return,DEM
endif

; initialize
tmin=0. & if ss eq 0 then tmin=10.^(min(logT))

; compute the normalization
if tmin eq 0 then norm=10.D^(alog10(EMtot)+alog10(ss)-ss*logTm) else $
norm = EMtot*ss/((10.D^(logTm))^(ss)-(tmin)^(ss))

; accumulate
for i=0L,ntmx-1L do begin
oo=where(tlog le logTm[i],moo)
if moo gt 0 then DEM[oo]=DEM[oo]+norm[i]*(10.D^(tlog[oo]))^(ss)
endfor

if vv gt 1000 then message,'HALT FOR DEBUGGING'

return,DEM
end