;+
; NAME:
;	IPS_WeightFnc
; PURPOSE:
;	Calculates weight function along the line of sight for IPS
; CATEGORY:
;	Physics
; CALLING SEQUENCE:
;	W = IPS_WeightFnc(ZAu, params)
; INPUTS:
;	ZAu			array[N]; type float
;					distance from Earth along line of sight (AU)
;	params		array[3]; type: float
;					params[0]: observing freq in MHz
;					params[1]: source size (arcsec)
;					params[2]: extra scale factor for source size
; OUTPUTS:
;	W			array[N]; type: float
;					weight factors
; CALLS:
;	IPS_params, nrSimpson
; PROCEDURE:
;	CALCULATE INTEGRAL FOR WEIGHTING FUNCTION OF Z.
;	SIMPSON METHOD DOUBLE INTEGRAL
;	SpectralFnc = q*fresnel(q)*size(q)*spectr(q)
;	fresnel= fresnel filter term
;	size   = source size term
;	spectr = spectrum fluctuation term
; MODIFICATION HISTORY:
;	MAR-1999, Paul Hick (pphick@ucsd.edu)
;		Adapted from STELab Fortran code
;-
FUNCTION SpectralFnc, q, arg

@compile_opt.pro		; On error, return to caller

THETA = arg[0]
WAVE  = arg[1]
Z     = arg[2]

Fnc = 0.0
IF q EQ 0.0 THEN RETURN, Fnc

;-------	Effect of source size
;-------	Spectrum fluctuation ~ q^-3

; Llimit of SpectralFnc integral by Q is determined by an exponent
; function of source size: 0.5*(q*THETA*Z)^2 = 10

F = q*THETA*Z
F = 0.5*F*F
IF F GE 10 THEN RETURN, Fnc

POWER = 3.0

Fnc = q*exp(-F)/q^POWER

F  = (q*q)*WAVE*Z/(4.0*!pi)	; Effect of fresnel filtering
DF = 2000.0*!pi
WHILE F GE 0.7 DO F -= DF
F = sin(F)

Fnc = F*F*Fnc

RETURN, Fnc  &  END


FUNCTION IPS_WeightFnc, ZAU, params

@compile_opt.pro		; On error, return to caller

IF n_elements(params) NE 3 THEN params = IPS_params(/nagoya)

FreqMHz		= params[0]
ArcSec		= params[1]
SourceScale	= params[2]

WAVE  = float( !physics.speedoflight )/(10*FreqMHz)	; Wavelength (km)
THETA = SourceScale*Arcsec/3600./!radeg				; Source size (Radians)
ZKM   = ZAU*float( !sun.au )*1.0E8					; LOS distance (km)

; limit of SpectralFnc integral by Q is determined by an exponent
; function of source size: 0.5*(q*THETA*Z)^2=10

WW = 0.0*ZAU
Xb = THETA*ZKM
UP = sqrt(60.0)/Xb
FOR I=0,n_elements(ZAU)-1 DO IF Xb[I] gt 0.0 THEN	$
	WW[I] = nrSimpson('SpectralFnc',0.,UP[I],eps=1.0E-5,maxn=13,arg=[THETA,WAVE,ZKM[I]])

RETURN, WW  &  END