%% whiscp.m calculates the numerically the integral used to estimate the 
%% equatorial electron density for whistlers.

function [op,int,dflm]=whiscp(oxdf,hedf,hydf,lsh,fn,nop)

% oxdf = fractional oxygen density at 1000km altitude (0->1)
% hedf = fractional helium density at 1000km altitude (0->1) 
% hydf = fractional hydrogen density at 1000km altitude (hydf=1-oxdf-hedf)
% lsh = L-shell of interest
% fn = 'nose frequency' (Hz)
% nop = number of points to break field line into
% op = integral output
% int = intergral terms
% dflm = distance along field lines

re=6371;                 %Earth radi (km)
r0=re+1000;              %Base radial distance (km) 
lonlsh=0;                %magnetic longitude to use (exact value unimportant)
k=1.38e-23;              %Boltzman constant (J/K) 
amu=1.66e-27;            %atomic mass unit (kg)

%Determine required co-ordinates
[lonmat,lonmatc,latmat,latmatc,latmatm,dismat,dflm]=flpar(lsh,lonlsh,r0,nop);

bf=magdip(dismat,latmatm);         %magnetic field at points on field line (T)
fh=28*1e9*bf;                      %electron gyro-frequency (Hz)
tmp=1600;                          %assumed temperature (K)
gv=9.8*(re/dismat(1))^2;           %gravity at r0 (m/s^2)
shhy=1e-3*k*tmp/(amu*gv);          %H+ scale height at r0 (km)
shhe=shhy/4;                       %He+ scale height at r0 (km)
shox=shhy/16;                      %O+ scale height at r0 (km)
we=2*pi/(24*3600);                 %angular velocity (rad/s)
z=dismat(1)-dismat(1)^2./dismat-(we^2/(2*gv))*((dismat.*cos(latmatm)).^2-(dismat(1)*cos(latmatm(1)))^2);

%Determine the I(s) parameter at each point along the field line
oxden=oxdf.*exp(-z./shox);             %O+ density variation
heden=hedf.*exp(-z./shhe);             %He+ density variation
hyden=hydf.*exp(-z./shhy);             %H+ density variation

neq=oxden(nop)+heden(nop)+hyden(nop);      %equatorial term
nf=oxden+heden+hyden;                      %density variation along field line
ist=sqrt(nf./neq);                         %I(s) term

int=sqrt(ist./(fn.*fh.*(1-fn./fh).^3));                          %integral term
op=sum(0.5*(int(1:(nop-1))+int(2:nop)).*diff(1e3*dflm));         %numerical sum of integral term

