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