%% polarfit calculates the coefficents of the quadratic used for the
%% polar profile. The profile is obtained by taking the Carpenter and
%% Anderson 1992 JGR midnight plasmatrough profile between L=3 and L=8
%% and calculating the distance along the field line from the base to
%% the equator of each point. 

clear all
close all

re=6371;
lsh=linspace(3,8,25);
sv=ones(size(lsh));
r0=re+350;

dflm=distfl2(0*sv,lsh,re*sv);
edenpt=1e6*5800*lsh.^(-4.5)+(1-exp(-(lsh-2)/10));
edenv=modvance(dflm);

c=polyfit(log10(dflm),log10(edenpt),2);

disp(' ')
fprintf(['c=[%6.3e %6.3e %6.3e];\n'],c)
disp(' ')

dflmf=linspace(3*re,10*re);
edf=10.^polyval(c,log10(dflmf));

dflm2=linspace(500,30e3);
edenv2=modvance(dflm2+re);
edenp=modper(dflm2+(re+350));
edf2=10.^polyval(c,log10(dflm2+re));
edeng=1e6*10.^(-3.09*log10(dflm2)+13.5);
lsh2=(dflm2+re)/re;
edenca=1e6*5800*lsh2.^(-4.5)+(1-exp(-(lsh2-2)/10));

figure(1)
h1=loglog(dflm,edenpt,'o',dflmf,edf);
grid
xlabel('Distance along field line (Re)')
ylabel('Log_{10} Plasmatrough Density 0 MLT')
legend('Carpenter and Anderson','Quadratic Fit',1)

figure(2)
h1=semilogx(edf2,dflm2/re+1,'o',edenv2,dflm2/re+1,'^',edenp,dflm2/re+1,'s',edeng,dflm2/re+1,'d',edenca,dflm2/re+1,'p');
grid
legend(h1,'Quadratic fit','Vance','Persoon','Gallagher','Carpenter',1)
