;+
; NAME : scattering_shift()
;
; PURPOSE : Estimates the radial shift, caused by radio-wave scattering, of the observed radio
; source away from its true location (and away from the Sun), using the simple,
; analytical method derived by Chrysaphi et al. (2018, ApJ, 868, 79).
;
; WRITTEN BY : Nicolina Chrysaphi @ Glasgow, UK
;
;
; CALLING SEQUENCE :
; result = scattering_shift(f=f, e2h=e2h [, /harmonic] [, /parker])
;
; INPUTS:
; f --- The observed frequency given in MHz. A 1D array of frequencies is also accepted.
; e2h --- The ratio of epsilon^2/h given in km^(-1).
; NOTE: The minimum and maximum values used by Chrysaphi et al. (2018, ApJ, 868, 79) for
; epsilon^2/h were 7e-5 and 4.5e-5 km^(-1), respectively.
;
; OPTIONAL INPUT (KEYWORD) PARAMETERS :
; harmonic --- If set, harmonic emission is assumed (i.e. f_observed = 2*f_pe).
; The default is harmonic=0, meaning that fundamental emission is assumed (i.e. f_observed = f_pe).
; parker --- If set, the modified Parker density model is used (see Kontar et al. (2019, ApJ, 884, 122) for details).
; The default is parker=0, meaning that the Newkirk (1961) coronal density model is assumed.
; NOTE: If f<25 MHz, use the Parker model. The Newkirk model does not suffice at larger distances.
;
; OUTPUT:
; Returns the radial shift of the observed radio source away from the true source, in solar radii.
; Note that scattering in the solar corona shifts sources away from the Sun (and the shift is frequency-dependent).
; The output has the same data type and dimensions as input f.
;
;
; OTHER NOTES :
; * This routine has also been translated into Python 3 and made publicly available.
; * Please acknowledge the use of this routine (and specify its version) in any publication, as well as
; cite Chrysaphi et al. (2018, ApJ, 868, 79) which describes the derivation of this analytical method.
; * The analytical method applied assumes that epsilon^2/h is constant with heliocentric distance, and
; takes the 1xNewkirk (1961) coronal density model (unless otherwise specified).
; * This method is valid for all radio emissions, irrespective of their exciter.
; * For a comprehensive explanation of the method and the assumptions made,
; see section 3.4 in Chrysaphi et al. (2018, ApJ, 868, 79), or
; see section 4.4 in Chrysaphi (2021), PhD thesis, University of Glasgow.
; * Access to the publication by Chrysaphi et al. (2018, ApJ, 868, 79):
; Link to the publisher's website: https://iopscience.iop.org/article/10.3847/1538-4357/aae9e5
; Link to NASA ADS : https://ui.adsabs.harvard.edu/abs/2018ApJ...868...79C
; Link to arXiv : https://arxiv.org/abs/1810.08026
; * Access to Chrysaphi (2021), PhD thesis, University of Glasgow:
; http://theses.gla.ac.uk/82010/
; * The scattering_shift() function calls the density_newkirk() and density_parker() functions which are included in this routine.
; * This routine was written in IDL (version 8.3.0).
;
; MODIFICATION HISTORY :
; v1.0 - 31-Oct-2019 - Nicolina Chrysaphi - Written as an IDL function.
; v1.1 - 14-Nov-2019 - Nicolina Chrysaphi - 1. Adjusted to accept an array of frequencies.
; 2. Added COMPILE_OPT IDL2 to all functions.
; v1.2 - 04-Jun-2020 - Nicolina Chrysaphi - Simplified by removing the R_newkirk() function.
; v2.0 - 25-Jan-2021 - Nicolina Chrysaphi - 1. Increased the range of distances calculated (up to 250 solar radii).
; 2. Added the option of using a modified Parker density model
; (as presented in Kontar et al. 2019, ApJ, 884, 122).
; v2.1 - 31-Mar-2021 - Nicolina Chrysaphi - 1. Added a user warning for the density model choice.
; 2. Routine was translated into Python 3 (and made publicly available).
;-
FUNCTION density_newkirk, r, N
;Newkirk (1961) coronal density model.
COMPILE_OPT IDL2
n_0 = 4.2e4
return, N*n_0*10.^(4.32/r)
;density in cm^-3
END
FUNCTION density_parker, r
;Modified Parker (1960) density model - see Kontar et al. (2019, ApJ, 884, 122) for details.
COMPILE_OPT IDL2
h1=20./960.
nc=3e11*exp(-(r-1.)/h1)
return, 4.8e9/r^14+ 3e8/r^6+1.39e6/r^2.3+nc
;density in cm^-3
END
FUNCTION scattering_shift, f=f, e2h=e2h, harmonic=harmonic, parker=parker
COMPILE_OPT IDL2
Nr = 29000L
r = findgen(Nr)/Nr*250+1.0
dr = r[1]-r[0]
IF keyword_set(parker) EQ 1 THEN BEGIN
dens = density_parker(r) ;cm^-3
ENDIF ELSE BEGIN
dens = density_newkirk(r, 1.0) ;cm^-3
ENDELSE
fpe = 8.9e-3*sqrt(dens) ;MHz
f = f ;MHz
eps2oh = e2h ;km^-1
Rsun = 6.96d5 ;km
FOR k=0,n_elements(f)-1 DO BEGIN
IF (f[k] LT 25.) AND (keyword_set(parker) NE 1) THEN BEGIN
print, 'WARNING: Choose your density model wisely...................'
stop
ENDIF
ENDFOR
scatt_shift = make_array(n_elements(f),/DOUBLE)
FOR i=0,n_elements(f)-1 DO BEGIN
fx = fpe[where(fpe LT f[i])]
rx = r[where(fpe LT f[i])]
IF keyword_set(harmonic) EQ 1 THEN BEGIN
tau_dr = (sqrt(!PI)/2.)*eps2oh*(fx^4/(( (2.*f[i])^2 - fx^2 )^2))
;For Harmonic emissions: f = 2*fx
ENDIF ELSE BEGIN
tau_dr = (sqrt(!PI)/2.)*eps2oh*(fx^4/((f[i]^2-fx^2)^2))
;For Fundamental emissions: f = fx
ENDELSE
tau = fltarr(n_elements(fx))
FOR j=0,n_elements(fx)-1 DO BEGIN
tau[j] = total(tau_dr[j:n_elements(rx)-1])*dr*Rsun
tau[n_elements(rx)-1] = tau[n_elements(rx)-2]
ENDFOR
tau1 = value_locate(tau, 1.0)
r1 = rx[tau1]
;Interpolate between the two closest tau=1 values to find exactly where tau=1:
INTEP, tau, rx, 1.0, r1_exact
scatt_shift[i] = ABS(r1_exact - min(rx))
;print, 'The radial shift away from the true source (of frequency ', f[i],' MHz) due to scattering effects was estimated to be = ', scatt_shift[i], ' R_solar .............'
ENDFOR
return, scatt_shift
END