;+ ; NAME: ; ThomsonLOSStep ; PURPOSE: ; Calculates integrated line-of-sight intensity for Thomson scattering. ; Integration is implemented as a sum over steps of equal size. ; CATEGORY: ; Physics: Thomson scattering ; CALLING SEQUENCE: FUNCTION ThomsonLOSStep, Pos_, Dir_, P, $ lower = Lower , $ upper = Upper , $ nstep = nstep , $ density = density , $ degrees = degrees , $ au = au , $ apm = APM , $ udark = udark , $ s10 = S10 ; INPUTS: ; Pos_ array[3,n]; type: float ; heliocentric location of observer: ; (ecliptic) longitude and latitude; and heliocentric ; distance ; (the 2nd dim can be absent, or can represent multiple ; dimensions) ; Dir_ array[2,m]; type: float ; topocentric (ecliptic) longitude and latitude of line of sight ; relative to Sun-observer direction ; (the 2nd dim can be absent, or can represent multiple ; dimensions) ; OPTIONAL INPUT PARAMETERS: ; /s10 if set intensities are returned in S10 units ; In this case the apparent magnitude APM MUST BE SPECIFIED ; udark=udark scalar; type: float; default: 0 ; limb darkening constant ; apm=APM scalar; type: float ; apparent magnitude of Sun ; density=density ; scalar; type: string; default: undefined ; name of function used to calculate the electron density. ; Keyword is passed to href=ThomsonSetupLOS=. ; lower=Lower ; scalar; type: float; default: 0 ; Lower limit of los integration (solar radii) ; If Lower <=0 then Lower = 0 is used. ; upper=Upper ; scalar; type: float; default: 9000 solar radii ; Upper limit of los integration (solar radii) ; If Lower <=0 then Upper = 9000.0 solar radii (~40 AU) is used. ; nstep=nstep ; scalar; type: integer; default: 100 ; number of elements in which to divide the range [Lower, Upper] ; ; /degrees if set all angles are in degrees (default: radians) ; /au if set all distances are in AU (default: solar radii) ; ; OUTPUTS: ; Result array[n,m]; type: double ; Integrated Thomson scattering intensity; ; units depend on setting of /S10: ; /S10 NOT set: per sterad in units of 10^-16 times the flux received ; from the solar disk (at the observer location ; /S10 set: S10 units ; ; P array[n,m]; type: float ; Polarization ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; InitVar, IsType, ToSolarRadii, ToRadians, gridgen, sphere_distance ; Distance2Sun, ThomsonBase, ThomsonSetupLOS, ThomsonLOSDensity ; ThomsonSolarFlux, ThomsonS10, destroyvar ; SEE ALSO: ; ThomsonLOSFar, ThomsonLOSRomb ; PROCEDURE: ; Densities are calculated by ThomsonLOSDensity. ; MODIFICATION HISTORY: ; SEP-1999, Paul Hick (UCSD/CASS) ; AUG-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ; Observer location can now also be an array. ; In addition the observer does not have to be in the ecliptic ; anymore. ;- InitVar, udark, 0 InitVar, S10, /key IF S10 THEN $ IF IsType(APM, /undefined) THEN $ message, 'no apparent magnitude for the Sun specified' r0 = 1.0d0/!sun.RAu rpau = ToSolarRadii(au=au) ; Converts to solar radii rpm = ToRadians(degrees=degrees);Converts to radians Tiny = 1.0d-5 Inf = 9000.0d0/rpau ; Maximum integration distance along line of sight InitVar, Lower, 0.0d0/rpau InitVar, Upper, Inf InitVar, nstep, 100 ; Check the integration limits; convert to solar radii L = Upper LT 0 S1 = (Lower > 0)*rpau S2 = ((L*Inf+(1B-L)*Upper) < Inf)*rpau ; nstep centers of integration intervals S = gridgen(nstep, range=[S1,S2], /open) ; Store the line of sight information into common block ; Get electron densities at all positions on all lines of sight pos = Pos_ los = Dir_ IF n_elements(pos) EQ 1 THEN pos = [0.0,0.0,pos] IF n_elements(dir) EQ 1 THEN los = [dir,0.0] szpos = size(pos) szlos = size(los) npos = szpos[szpos[0]+2]/szpos[1] nlos = szlos[szlos[0]+2]/szlos[1] pos = reform(pos,szpos[1],npos,/overwrite) ; [3,npos] los = reform(los,szlos[1],nlos,/overwrite) ; [2,nlos] pos[0:1,*] *= rpm ; Radians pos[ 2,*] *= rpau ; Solar radii los *= rpm ; Radians ThomsonSetupLOS, pos, los, density=density Den = ThomsonLOSDensity(S) ; # electrons/cm^3 pos = SuperArray(pos,nlos,after=1) ; [3,nlos,npos] pos = reform(pos,szpos[1],nlos*npos,/overwrite); [3,nlos*npos] los = SuperArray(los,npos,/trail) ; [2,nlos,npos] los = reform(los,szlos[1],nlos*npos,/overwrite); [2,nlos*npos] ; Calculate intensity per electron at all ; positions on all lines of sight ; Direction to sun (in same sun-centered coordinate system as los, ; i.e. the sun is at lng=0 at the opposite latitude) tmp = pos[0:1,*] ; [2,nlos*npos] tmp[0,*] = 0.0d0 tmp[1,*] *= -1.0d0 ; Elongations for all lines of sight elo = sphere_distance(tmp,los) ; [nlos*npos] pos = SuperArray(pos,nstep,after=1) ; [3,nstep,nlos*npos] elo = SuperArray(elo,nstep,/lead) ; [ nstep,nlos*npos] pos = reform(pos,szpos[1],nstep*nlos*npos,/overwrite); [3,nstep*nlos*npos] elo = reform(elo, nstep*nlos*npos,/overwrite); [ nstep*nlos*npos] ; Get distance to Sun and SinChi angle rsun = reform(pos[2,*]) ; [nstep*nlos*npos] tmp = SuperArray(S,nlos*npos,/trail) tmp = reform(tmp,nstep*nlos*npos,/overwrite); [nstep*nlos*npos] tmp = rsun*Distance2Sun(elo, tmp/rsun, SinChi) tmp = ThomsonBase(tmp,SinChi,udark,P,It,Itr) ; Sum over the first dimension (along line of sight) ; Save the structure of the input array Dir with LOS directions It = total(Den*It ,1) Itr = total(Den*Itr,1) It = It+It-Itr ; Total intensity P = Itr/It n = nstep-1 dS = (S[n]-S[0])/n ; Step size (solar radii) cv = 0.5d0*!physics.thomsoncrosssection ; Replacing RSun by r0 if S10=1 doesn't change the result, but ; is faster when RSun is an array. rsun = reform(rsun,nstep,nlos,npos,/overwrite) rsun = reform(rsun[0,*,*]) CASE S10 OF 0: cv /= ThomsonSolarFlux(rsun,udark) ; Intensity/Solar Flux 1: cv *= 1.0d-16*ThomsonS10(r0,APM)/ThomsonSolarFlux(r0 ,udark) ; S10 units ENDCASE cv *= !sun.R*dS*It destroyvar, tmp IF szlos[0] GT 1 THEN boost, tmp, szlos[2:szlos[0]] IF szpos[0] GT 1 THEN boost, tmp, szpos[2:szpos[0]] CASE IsType(tmp,/defined) OF 0: BEGIN cv = cv[0] P = P [0] END 1: BEGIN cv = reform(cv,tmp,/overwrite) P = reform(P ,tmp,/overwrite) END ENDCASE RETURN, cv & END