;+ ; NAME: ; IntegrateLOS ; PURPOSE: ; Calculate a sky map by integrating through a 3D heliospheric matrix. ; Two forms of integrated expressions are implemented: ; (1) A weighted integration (used for g-level and IPS velocity): ; F = Integral(F*W ds)/Integral(W ds) ; (2) Straight integration (used for Thomson scattering): ; F = Integral(F ds) ; See PROCEDURE for more details. ; CATEGORY: ; sat/idl/util ; CALLING SEQUENCE: FUNCTION IntegrateLOS, R, FF, WW, $ ; Keywords passed to InterpolateHeliosphere R3D = R3D , $ dR3D = dR3D , $ xcrange = xcrange , $ cv2carrington=cv2carrington,$ xlrange = xlrange , $ opengrid= opengrid , $ xcgrid = xcgrid , $ xlgrid = xlgrid , $ rrgrid = rrgrid , $ fillbad = fillbad , $ degrees = degrees , $ minRR = minRR , $ f3dfnc = f3dfnc , $ f3darg = f3darg , $ w3dfnc = w3dfnc , $ w3darg = w3darg , $ rr_earth= rr_earth , $ ut_earth= ut_earth , $ pa_earth= pa_earth , $ elo_earth=elo_earth , $ elo_sun = elo_sun , $ silent = silent ; INPUTS: ; R array[3,*,n]; type: float ; Locations in spherical coordinates for all points along all ; lines of sight. '*' can be absent or represent more than one dimension. ; E.g. for a typical skymap there will be two dimensions l,m ; representing longitude and latitude. 'n' represents the ; 'depth' dimension along the lines of sight. ; ; R[0,*,n] contains longitude angles in [0,360] ; R[1,*,n] contains latitudes in [-90,90] ; R[2,*,n] contains heliocentric distances in AU ; ; The last dimension is the 'depth' dimension along the line of sight over ; which the integration is performed. ; ; 'Densities' and 'Weights' can be specified in two ways. ; ; If R3D and xcgrid are NOT specified then FF and WW directly specify densities ; and weights along all lines of sight: ; ; FF array[*,n]; type: float ; volume data at locations 'R'. ; WW array[*,n]; type: float ; weights at locations 'R'. ; ; If R3D or xcgrid are specified than calls to InterpolateHeliosphere ; are made with arrays FF and WW to determine the function values and ; weight along each line of sight. ; ; In this case the following keywords are passed to InterpolateHeliosphere ; R ; FF or WW ; R3D = R3D ; dR3D = dR3D ; xcrange = xcrange ; cv2carrington=cv2carrington ; xlrange = xlrange ; opengrid= opengrid ; xcgrid = xcgrid ; xlgrid = xlgrid ; rrgrid = rrgrid ; fillbad = fillbad ; ; See href=InterpolateHeliosphere for more information. ; NOTE: R[0,*] MUST BE longitude angles, even though InterpolateHeliosphere ; in principle allows both R[0,*] and xcrange to be Carrington variables. ; R[0,*] must be a longitude angle because R is also passed as ; argument to f3dfnc and w3dfnc, where it is used in combination ; with the longitude angles in rr_earth ; ; OPTIONAL INPUT PARAMETERS: ; /degrees if set then all angles are in degrees (default: radians) ; ; The integrands are set up using a number of additional keywords: ; ; ut_earth=ut_earth ; rr_earth=rr_earth ; position Earth in same coordinates as R. ; pa_earth=pa_earth ; array[*,n]; type: float ; line of sight position angle measured counterclockwise ; from ecliptic north ; elo_earth=elo_earth ; array[*,n]; type: float ; angle Sun-observer-segment on line of sight (i.e. the elongation) ; for all line of sight segments (note that array[*,0] = array[*,1] etc. ; elo_sun=elo_sun ; array[*,n]; type: float ; angle observer-Sun-segment on line of sight ; OUTPUTS: ; Result array[*]; type: float ; sky map with line-of-sight integrated quantities ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; IsType, InitVar, InterpolateHeliosphere, BadValue ; PROCEDURE: ; > Two types of expressions are evaluated ; ; (1) = Integral[ w3dfnc(s,W(s))*f3dfnc(s,F(s)) ds]/ Integral[ w3dfnc(s,W(s)) ds] ; and ; (2) = Integral[ w3dfnc(s,W(s))*f3dfnc(s,F(s)) ds] ; ; The line of sight integrations is done for all lines of sight specified. ; The integration dimension always is the last dimension in the input arrays. ; The integration is performed simply by summing over the integration dimension, ; implicitly assuming a constant stepsize ds: ; ; (1) = Sum[ w3dfnc(s,W(s))*f3dfnc(s,F(s)) ]/ Sum[ w3dfnc(s,W(s)) ] ; and ; (2) = Sum[ w3dfnc(s,W(s))*f3dfnc(s,F(s)) ] ; ; Note that for expression (1) the stepsize ds will cancel. For expression (2) ; it is the users responsibility to absorb the stepsize ds into one of the ; input arrays F or W. ; ; > W(s) and F(s) represent two physical quantities known on a heliographic grid ; (for example solar wind velocity and density, and the s-dependence reflects some ; geometrical aspect of the integration. ; ; In all applications I can think of the dependence of f3dfnc on F(s) will be linear: ; f3dfnc(s,t) = f3dfnc(s)*t. In that case the integration looks like: ; = Integral[ w3dfnc(s,W(s))*f3dfnc(s)*F(s) ds ]/ Integral[ w3dfnc(s,W(s)) ds ] ; This form applies to the calculation of IPS velocities, where the weight factor ; depends on the density: w3dfnc is 'href=IPS_velocity_w=', W(s) = n(s), and ; the perpendicular velocity is integrated: f3dfnc(s) = sin(angle(s)) (is ; 'href=IPS_velocity=') and F(s) = V(s). ; ; In its simplest form w3dfnc(s,t) = w3dfnc(s) and f3dfnc(s,t) = t. In that case: ; = Integral[ w3dfnc(s)*F(s) ds ]/ Integral[ w3dfnc(s) ds ] ; i.e. is a weighted average of F(s) along the line of sight. This form applies ; to the calculation of IPS g-level with a weight factor that depends only on the ; line of sight geometry: w3dfnc is 'href=IPS_hlevel_w=', and the quantity F(s) is ; the 'gamma function' (introduced in the tomography manifesto). ; ; A special case occurs when there is no denominator: ; = Integral[ w3dfnc(s)*F(s) ds ] ; This applies to the Thomson scattering brightness where w3dfnc is the scattered ; intensity from a single electron ('href=ThomsonBrightness='), and F(s) is the ; solar wind density n(s). The stepsize ds is absorbed into ThomsonBrightness. ; ; > Bad values in the input array F3D are detected with the IDL 'finite' function. ; This is done in InterpolateHeliosphere. The treatment of these bad values ; depends on the setting of 'fillbad'. ; ; MODIFICATION HISTORY: ; AUG-1999, Paul Hick (UCSD/CASS) ; AUG-2004, Paul Hick (UCSD/CASS) ; Bug fix: w3dfnc was not be used if WW did not exist. This is wrong ; when w3dfnc only uses R and not WW (as for the g-level determination). ; Bug fix: if neither w3dfnc nor WW is defined (as for Thomson ; scattering the average rather than the total along the line of ; sight was calculated. ; Also rearranged some of the arrangement. Now the presence of R3D ; triggers a call to InterpolateHeliosphere. ; JUN-2006, Paul Hick (UCSD/CASS) ; Fixed bug in the summation expressions. It was explicitly set to ; the third dimension. Now it is set to the last dimension present. ; in F and W. Added grid keywords xcgrid, xlgrid and rrgrid to call to ; InterpolateHeliosphere. ; FEB-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ; Added test for existence of xcgrid to trigger call to ; InterpolateHeliosphere. ;- InitVar, silent, 0 sample = IsType(R3D,/defined) OR IsType(xcgrid,/defined) ; Get the density values along all lines of sight, either by interpolation ; (sample SET) or directly from input argument FF (sample NOT set). CASE IsType(FF, /defined) OF 0: IF silent LE 0 THEN message, /info, 'no density values specified' 1: F = sample ? InterpolateHeliosphere(R, FF, R3D, dR3D, $ xcrange=xcrange, xlrange=xlrange, opengrid=opengrid, fillbad=fillbad, $ xcgrid=xcgrid, xlgrid=xlgrid, rrgrid=rrgrid, $ cv2carrington=cv2carrington, degrees=degrees) : FF ENDCASE ; Apply f3dfnc to densities F. Note that technically F may not exist (and f3dfnc only ; would operate on positional info rr_earth and F). IF IsType(f3dfnc,/defined) THEN $ F = call_function(f3dfnc, rr_earth, R, F, f3darg, degrees=degrees, $ ut_earth=ut_earth, pa_earth=pa_earth, elo_earth=elo_earth, elo_sun=elo_sun) CASE IsType(WW, /defined) OF 0: IF silent LE 0 THEN message, /info, 'no weight values specified' 1: W = sample ? InterpolateHeliosphere(R, WW, R3D, dR3D, $ xcrange=xcrange, xlrange=xlrange, opengrid=opengrid, fillbad=fillbad, $ xcgrid=xcgrid, xlgrid=xlgrid, rrgrid=rrgrid, $ cv2carrington=cv2carrington, degrees=degrees) : WW ENDCASE ; Apply w3dfnc to weights W. Note that technically W may not exist (and w3dfnc ; only would operate on positional info rr_earth and F). This actually happens ; for g-level calculations, when w3dfnc='IPS_hlevel_w' does not use W. IF IsType(w3dfnc,/defined) THEN $ W = call_function(w3dfnc, rr_earth, R, W, w3darg, degrees=degrees, $ ut_earth=ut_earth, pa_earth=pa_earth, elo_earth=elo_earth, elo_sun=elo_sun) ; If W is not defined at this point, it was not specified as input and no ; w3dfnc function was specified. In this case the line of sight densities ; in F are just summed. This is used to calculate the Thomson scattering ; brightness when FF is the solar wind density and f3dfnc='ThomsonBrightness'. CASE IsType(W,/defined) OF 0: BEGIN goodF = finite(F) bad = where(1-goodF) IF bad[0] NE -1 THEN BEGIN F[bad] = 0.0 goodF[bad] = 1B ENDIF F = total(F,size(F,/n_dim),/nan) ; Sum over last dimension END 1: BEGIN goodF = finite(F) goodW = finite(W) bad = where(goodF AND NOT goodW) IF bad[0] NE -1 THEN BEGIN F[bad] = BadValue(F) goodF[bad] = 0B ENDIF bad = where(goodW AND NOT goodF) IF bad[0] NE -1 THEN BEGIN F[bad] = 0.0 goodF[bad] = 1B ;W[bad] = BadValue(W) ;goodW[bad] = 0B ENDIF F = total(F*W,size(F,/n_dim),/nan)/total(W,size(W,/n_dim),/nan) END ENDCASE InitVar, minRR, 0 W = where(total(goodF,size(goodF,/n_dim)) LT minRR) IF W[0] NE -1 THEN F[W] = BadValue(F) RETURN, F & END