%% pedcon.m estimates the Pedersen conductivity along the specified
%% field lines in both hemisphere. The Pedson formula used is given
%% in Brekke eqn 7.55, p324 and for the ion-neutral and electron-
%% neutral collision frequencies uses the approximate equations given
%% by Brekke eqn 7.11 and 7.13. MSIS-90 and IRI95 are used to 
%% calaculate the required neutral and electron densities.
%% The integrated Pedersen conductivity in both hemispheres is returned
%% as the same sized array as input, others parameters are returned in
%% format where each column corresponds to a input value.
%% NOTE - should not be used below L=1.5.
%% NOTE (9/1/00) replace Brekke eqn 7.13 with electron neutral collsions for
%% N2, O2 and O (major species in this height range) from Table 3, Schunk 
%% and Nagy, 1980.

function [pedn,peds,altv,pedfn,pedfs,cfinn,cfins,cfenn,cfens]=pedcon(lsh,lonlsh,timmat,day,month,year)

% lsh = L-shell of interest 
% lonlsh = magnetic longitude of field line (rads)
% timmat = UT of interest (fractional hours)
% day = day of interest
% month = month of interest
% year = year of interest
%
% pedn = integrated northern hemisphere Pedersen conductivty (S)
% peds = integrated southern hemisphere Pedersen conductivty (S)
% altv = altitude vector used (km)
% pednf = northern hemisphere Pedersen conductivty values (S)
% pedfs = southern hemisphere Pedersen conductivty values (S)
% cfinn = northern hemisphere ion-neutral collision frequency (s^-1)
% cfins = southern hemisphere ion-neutral collision frequency (s^-1)
% cfenn = northern hemisphere electron-neutral collision frequency (s^-1)
% cfens = southern hemisphere electron-neutral collision frequency (s^-1)

re=6371;                                %Earth radi (km)
amu=1.66e-27;                           %atomic mass unit (kg)
em=9.1e-31;                             %electron mass (kg)
ec=1.6e-19;                             %electronic charge (C)
bo=30.4e-6;                             %equatorial magnetic field (T)

sz=size(lsh);                           %size of input

lsh=lsh(:);
lonlsh=lonlsh(:);
timmat=timmat(:);
day=day(:);
month=month(:);
year=year(:);

dno=daynum(day,month,year);             %day number

%Determine required solar parameters
[kp,kpmin,apday,f107,f107a]=solarval(timmat,dno,year);

%Calculate the Pederson values between 90km and 600km. Set up an array
%for each field line between these altitudes.
altv=[90:5:150 160:10:300 320:20:400];            %altitude vector (km)

nopy=length(lsh);
nopx=length(altv);
svx=ones(1,nopx);
svy=ones(nopy,1);

latmatm=lshpar(lsh*svx,-svy*svx,svy*(re+altv));   %required magnetic latitude values
[lonn,latn]=mag2geo3(lonlsh*svx,latmatm);         %northern geographic co-ordinates
[lons,lats]=mag2geo3(lonlsh*svx,-latmatm);        %southern geographic co-ordinates

%Calculate required ion densities (percentages)
[eden,hyden,oxden,heden,o2den,noden,tmpe]=modiri7([lonn;lons],[latn;lats],[timmat;timmat]*svx,[day;day]*svx,[month;month]*svx,[year;year]*svx,[svy*altv;svy*altv],0);

%Replace -1 electron temperature values with first non -1 value
[ti,tj]=find(tmpe==-1);
poi=find(tmpe==-1);
tmpe2=tmpe;
tmpe2(tmpe2==-1)=nan;                             %reset -1 to NaN
tmin=min(tmpe2,[],2);                             %minimum temperature for each field line
tmpe(poi)=tmin(ti);                               %replace required temperatures

%Calculate required neutral densities 
[hent,oxnt,n2nt,o2nt,arnt,hynt,nint]=modmsis7([lonn;lons],[latn;lats],[svy*altv;svy*altv],[timmat;timmat]*svx,[day;day]*svx,[month;month]*svx,[year;year]*svx,[f107a;f107a]*svx,[f107;f107]*svx,[apday;apday]*svx);

nnt=hent+oxnt+n2nt+o2nt+arnt+hynt+nint;                               %total neutral density (m^-3)
mnm=(4*hent+16*oxnt+28*n2nt+32*o2nt+18*arnt+hynt+14*nint)./nnt;       %mean neutral molecular mass (amu)

%Calculate the magnetic field at the points of interest
bf=bo*(re./(svy*(altv+re))).^3.*sqrt(1+3*sin(latmatm).^2);  %magnetic flux (T)

ni=1:nopy;                                                  %northern index terms
si=(nopy+1):(2*nopy);                                       %southern index terms

cfinn=2.6e-15*(nnt(ni,:)+eden(ni,:))./sqrt(mnm(ni,:));      %northern ion-neutral collision frequency (s^-1)
cfins=2.6e-15*(nnt(si,:)+eden(si,:))./sqrt(mnm(si,:));      %southern ion-neutral collision frequency (s^-1)

cfen2=2.33e-17*n2nt.*(1-1.21e-4*tmpe).*tmpe;                %N2 electron neutral collision frequency (s^-1)
cfeo2=1.82e-16*o2nt.*(1+3.6e-2*sqrt(tmpe)).*sqrt(tmpe);     %O2 electron neutral collision frequency (s^-1)
cfeox=8.9e-17*oxnt.*(1+5.7e-4*tmpe).*sqrt(tmpe);            %O electron neutral collision frequency (s^-1)

cfenn=cfen2(ni,:)+cfeo2(ni,:)+cfeox(ni,:);                  %northern electron-neutral collision frequency (s^-1)
cfens=cfen2(si,:)+cfeo2(si,:)+cfeox(si,:);                  %southern electron-neutral collision frequency (s^-1)

%cfenn=5.4e-16*nnt(ni,:).*sqrt(tmpe(ni,:));                  %northern electron-neutral collision frequency (s^-1)
%cfens=5.4e-16*nnt(si,:).*sqrt(tmpe(si,:));                  %southern electron-neutral collision frequency (s^-1)

gfin=ec.*bf./(amu*mnm(ni,:));                               %northern ion gyrofrequency (Hz)
gfis=ec.*bf./(amu*mnm(si,:));                               %southern ion gyrofrequency (Hz)
gfe=ec.*bf./em;                                             %electron gyrofrequency (Hz)

kin=gfin./cfinn;                                            %northern ion mobility co-efficent
kis=gfis./cfins;                                            %southern ion mobility co-efficent
ken=gfe./cfenn;                                             %northern electron mobility co-efficent
kes=gfe./cfens;                                             %southern electron mobility co-efficent

pedfn=(eden(ni,:)*ec./bf).*(ken./(1+ken.^2)+kin./(1+kin.^2));    %northern Pederson conductivity (S/m)
pedfs=(eden(si,:)*ec./bf).*(kes./(1+kes.^2)+kis./(1+kis.^2));    %southern Pederson conductivity (S/m)

altm=1e3*svy*diff(altv);
poi1=1:(nopx-1);
poi2=2:nopx;
pedn=sum(0.5*(pedfn(:,poi2)+pedfn(:,poi1)).*altm,2);        %northern height integrated Pederson conductivity (S)
peds=sum(0.5*(pedfs(:,poi2)+pedfs(:,poi1)).*altm,2);        %southern height integrated Pederson conductivity (S)

