%% plasieq.m calculates and plots the equatorial e- and ion densities,
%% as well as the TEC and mean ionospheric height in the equatorial 
%% plane using data files produced by plasimod1.m. 
%%
%% 10/11/2001 - update to include He+.
%% 5/3/2002 - remove different files name extension for IRI CCIR or URSI parameters (store iricf in arrays).
%% 5/3/2002 - He+ inclusion now optional. 
%%

clear

timmat=input('UT of interest? ');
year=input('Year of interest? ');
di=input('Day of interest (negative value for day number)? ');
if di<0
 [day,month]=daynum2(-di,year);
 dno=-di;
else 
 month=input('Month of interest? ');
 day=di;
 dno=daynum(di,month,year);
end 

%Determine the maximum allowed altitude.
re=6371;                                %Earth radi (km)
lshcp=plasires2(dno,year);
altmx=1e3*floor((lshcp-1)*re/1e3);

altmax=input(['Maximum altitude (max ',int2str(altmx),' km)? ']);
tecalt=input('Input vector with altitude band values required \n for TEC (0=just total TEC, 5 values max.)? ');

[eden,oxden,hyden,heden,lonmatm,dismat,mltv,hes]=plasieq1(altmax,timmat,dno,year);

kp=solarval(timmat,dno,year);
kp=round(10*kp)/10;
svy=ones(size(dismat,1),1);
mltm=(pi/12)*(svy*mltv);
sz=size(mltm);
mltm(:,sz(2))=0;
[xmat,ymat]=pol2cart(mltm,dismat);
re=6371;
tecalt=tecalt(:)';
disv=dismat(:,1);
[tec,mih,tectot]=teccal1(eden',disv',tecalt);
techy=teccal1(hyden',disv',0);
if hes==1
 teche=teccal1(heden',disv',0);
else
 teche=nan*ones(size(techy));
end 
[hytrh,hytrd]=tranht(dismat'-re,oxden',hyden');
if hes==1
 [hetrh,hetrd]=tranht(dismat'-re,oxden',heden');
else
 hetrh=nan*ones(size(hytrh));
 hetrd=nan*ones(size(hytrd));
end
tectr=teccal2(eden',dismat(:,1)',hytrh);

timmat2=round(10*timmat)/10;
dismx=max(disv)/re;

figure(1)
pcolor(xmat/re,ymat/re,log10(eden))
axis equal
shading flat
colorbar
hold on
[xme1,yme1]=pol2cart((pi/12)*linspace(-6,6,50),1);
[xme2,yme2]=pol2cart((pi/12)*linspace(6,18,50),1);
fill(xme1,yme1,'k')
fill(xme2,yme2,'g')
hold off
titl=['UT=',num2str(timmat2),'  ',int2str(day),'/'int2str(month),'/',int2str(year),' Kp=',num2str(kp)];
title(['Log_{10} e^-  ',titl],'fontsize',15,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold')
text(0.6*dismx,0,'0 MLT','fontsize',14,'fontweight','bold','color','w')
text(0,0.8*dismx,'6 MLT','fontsize',14,'fontweight','bold','color','w')
text(-0.8*dismx,0,'12 MLT','fontsize',14,'fontweight','bold','color','w')
text(0,-0.8*dismx,'18 MLT','fontsize',14,'fontweight','bold','color','w')
xlabel('Earth radi (L-shell)','fontsize',12,'fontweight','bold')
ylabel('Earth radi (L-shell)','fontsize',12,'fontweight','bold')

oxden2=oxden;
oxden2(oxden2==0)=NaN;
hyden2=hyden;
hyden2(hyden2==0)=NaN;
heden2=heden;
heden2(heden2==0)=NaN;

figure(2)
pcolor(xmat/re,ymat/re,log10(oxden2))
axis equal
shading flat
colorbar
hold on
fill(xme1,yme1,'k')
fill(xme2,yme2,'g')
hold off
title(['Log_{10} 0^+  ',titl],'fontsize',15,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold')
text(0.6*dismx,0,'0 MLT','fontsize',14,'fontweight','bold','color','w')
text(0,0.8*dismx,'6 MLT','fontsize',14,'fontweight','bold','color','w')
text(-0.8*dismx,0,'12 MLT','fontsize',14,'fontweight','bold','color','w')
text(0,-0.8*dismx,'18 MLT','fontsize',14,'fontweight','bold','color','w')
xlabel('Earth radi (L-shell)','fontsize',12,'fontweight','bold')
ylabel('Earth radi (L-shell)','fontsize',12,'fontweight','bold')

figure(3)
pcolor(xmat/re,ymat/re,log10(hyden2))
axis equal
shading flat
colorbar
hold on
fill(xme1,yme1,'k')
fill(xme2,yme2,'g')
hold off
title(['Log_{10} H^+  ',titl],'fontsize',15,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold')
text(0.6*dismx,0,'0 MLT','fontsize',14,'fontweight','bold','color','w')
text(0,0.8*dismx,'6 MLT','fontsize',14,'fontweight','bold','color','w')
text(-0.8*dismx,0,'12 MLT','fontsize',14,'fontweight','bold','color','w')
text(0,-0.8*dismx,'18 MLT','fontsize',14,'fontweight','bold','color','w')
xlabel('Earth radi (L-shell)','fontsize',12,'fontweight','bold')
ylabel('Earth radi (L-shell)','fontsize',12,'fontweight','bold')

if hes==1
 figure(4)
 pcolor(xmat/re,ymat/re,log10(heden2))
 axis equal
 shading flat
 colorbar
 hold on
 fill(xme1,yme1,'k')
 fill(xme2,yme2,'g')
 hold off
 title(['Log_{10} He^+  ',titl],'fontsize',15,'fontweight','bold')
 set(gca,'fontsize',12,'fontweight','bold')
 text(0.6*dismx,0,'0 MLT','fontsize',14,'fontweight','bold','color','w')
 text(0,0.8*dismx,'6 MLT','fontsize',14,'fontweight','bold','color','w')
 text(-0.8*dismx,0,'12 MLT','fontsize',14,'fontweight','bold','color','w')
 text(0,-0.8*dismx,'18 MLT','fontsize',14,'fontweight','bold','color','w')
 xlabel('Earth radi (L-shell)','fontsize',12,'fontweight','bold')
 ylabel('Earth radi (L-shell)','fontsize',12,'fontweight','bold')
end


mltv2=mltv(:);
cv=['r';'g';'b';'c';'m';'y'];
cv2=cv(1:(length(tecalt)+1));
cv2=flipud(cv2);
for i=1:6
 eval(['legt',int2str(i),'=[];'])
end 
if tecalt==0
 legt1='Total TEC';
else 
 tecalt2=[0 tecalt altmax]; 
 for j=1:(length(tecalt)+1)
  legt0=[num2str(tecalt2(j)),' to ',num2str(tecalt2(j+1)),' km'];
  eval(['legt',int2str(j),'=legt0;']);
 end 
end
figure(5)
h0=area(mltv2,tec);
if tecalt==0
 set(h0,'facecolor',cv(1))
else
 for j=1:(length(tecalt)+1)
  set(h0(j),'facecolor',cv(j))
 end 
end 
grid
axis([0 24 0 10*ceil(max(tectot)/10)])
xlabel('Magnetic Local Time','fontsize',12,'fontweight','bold')
ylabel('Total Electron Content (TECU)','fontsize',12,'fontweight','bold')
title(['TEC Altitude bands: ',titl],'fontsize',15,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold','layer','top')
if ~isempty(legt6)
 legend(h0,legt1,legt2,legt3,legt4,legt5,legt6)
elseif ~isempty(legt5)
 legend(h0,legt1,legt2,legt3,legt4,legt5)
elseif ~isempty(legt4)
 legend(h0,legt1,legt2,legt3,legt4)
elseif ~isempty(legt3)
 legend(h0,legt1,legt2,legt3)
elseif ~isempty(legt2)
 legend(h0,legt1,legt2)
else
 legend(h0,legt1)
end 

for i=1:6
 eval(['legt',int2str(i),'=[];'])
end 
if tecalt==0
 legt1='Total TEC';
else 
 tecalt2=[altmax fliplr(tecalt) 0];
 for j=1:(length(tecalt)+1)
  legt0=[num2str(tecalt2(j+1)),' to ',num2str(tecalt2(j)),' km'];
  eval(['legt',int2str(j),'=legt0;']);
 end 
end
figure(6)
h0=area(mltv2,fliplr(tec));
if tecalt==0
 set(h0,'facecolor',cv2(1))
else
 for j=1:(length(tecalt)+1)
  set(h0(j),'facecolor',cv2(j))
 end 
end 
grid
axis([0 24 0 10*ceil(max(tectot)/10)])
xlabel('Magnetic Local Time','fontsize',12,'fontweight','bold')
ylabel('Total Electron Content (TECU)','fontsize',12,'fontweight','bold')
title(['TEC Altitude bands: ',titl],'fontsize',15,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold','layer','top')
if ~isempty(legt6)
 legend(h0,legt1,legt2,legt3,legt4,legt5,legt6)
elseif ~isempty(legt5)
 legend(h0,legt1,legt2,legt3,legt4,legt5)
elseif ~isempty(legt4)
 legend(h0,legt1,legt2,legt3,legt4)
elseif ~isempty(legt3)
 legend(h0,legt1,legt2,legt3)
elseif ~isempty(legt2)
 legend(h0,legt1,legt2)
else
 legend(h0,legt1)
end 

figure(7)
hp=plot(mltv2,mih,'bo-');
grid
xlabel('Magnetic Local Time','fontsize',12,'fontweight','bold')
ylabel('Mean ionospheric height (km)','fontsize',12,'fontweight','bold')
title(['Mean Ionospheric Height: ',titl],'fontsize',15,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold')
set(hp,'linewidth',2)
set(hp(1),'markerfacecolor','b')
axis([0 24 0 1.1*round(max(mih))])

figure(8)
subplot(1,2,1)
hp=plot(mltv2,hytrh,'bo-',mltv2,hetrh,'r^-');
grid
title('Transition Heights','fontsize',12,'fontweight','bold')
xlabel('Magnetic Local time','fontsize',10,'fontweight','bold')
ylabel('Transition Height (km)','fontsize',10,'fontweight','bold')
set(gca,'fontsize',10,'fontweight','bold')
set(hp,'linewidth',2)
set(hp(1),'markerfacecolor','b')
set(hp(2),'markerfacecolor','r')
axis([0 24 0 1.1*round(max([hytrh;hetrh]))])
legend(hp,'O^+/H^+','O^+/He^+',1)  
subplot(1,2,2)
hytrd(hytrd==0)=NaN;
hetrd(hetrd==0)=NaN;
hp=plot(mltv2,log10(hytrd),'bo-',mltv2,log10(hetrd),'^r-');
grid
xlabel('Magnetic Local Time','fontsize',10,'fontweight','bold')
ylabel('Transition Density(log10)','fontsize',10,'fontweight','bold')
set(gca,'fontsize',10,'fontweight','bold')
set(hp,'linewidth',2)
set(hp(1),'markerfacecolor','b')
set(hp(2),'markerfacecolor','r')
axis([0 24 0.975*log10(round(min([hytrd;hetrd]))) 1.025*log10(round(max([hytrd;hytrd])))])
legend(hp,'O^+/H^+','O^+/He^+',1)  

figure(9)
subplot(2,1,1)
h0=area(mltv2,tectr);
grid
axis([0 24 0 10*ceil(max(tectot)/10)])
xlabel('Magnetic Local Time','fontsize',12,'fontweight','bold')
ylabel('Total Electron Content (TECU)','fontsize',12,'fontweight','bold')
title(['TEC above and below O^+/H^+ transition height: '],'fontsize',15,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold','layer','top')
legend(h0,'TEC below','TEC above')
subplot(2,1,2)
h0=plot(mltv2,tectr(:,2),'r^-');
grid
axis([0 24 0 10*ceil(max(tectr(:,2))/10)])
xlabel('Magnetic Local Time','fontsize',12,'fontweight','bold')
ylabel('Total Electron Content (TECU)','fontsize',12,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold')
set(h0,'linewidth',2,'markerfacecolor','r')
legend(h0,'TEC above')

figure(10)
h0=plot(mltv2,100*techy./tectot,'bo-',mltv2,100*teche./tectot,'r^-');
grid
axis([0 24 0 10*ceil(max(100*techy./tectot)/10)])
title('Percentage of TEC due to H^+ and He^+','fontsize',15,'fontweight','bold')
xlabel('Magnetic Local Time','fontsize',12,'fontweight','bold')
ylabel('H^+ % of TEC','fontsize',12,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold')
set(h0,'linewidth',2)
set(h0(1),'markerfacecolor','b')
set(h0(2),'markerfacecolor','r')
legend(h0,'H^+','He^+',1)

lshm=dismat/re;
[eca,eptca]=modca(lshm(:,1),year*svy,month*svy,day*svy,timmat*svy);
edeng=modgall2(lshm(:,1),year*svy,month*svy,day*svy,timmat*svy);
eca(eca==0)=nan;
eptca(eptca==0)=nan;
edeng(edeng==0)=nan;
edi1=find(mltv==2);
edi2=find(mltv==8);
edi3=find(mltv==14);
edi4=find(mltv==20);
figure(11)
hp=semilogy(lshm(:,1),eden(:,edi1),'gs-',lshm(:,1),eden(:,edi2),'yo-',lshm(:,1),eden(:,edi3),'mp-',lshm(:,1),eden(:,edi4),'c^-',lshm(:,1),1e6*eca,'b',lshm(:,1),1e6*eptca,'r',lshm(:,1),1e6*edeng,'k');
set(hp(1),'markerfacecolor','g','markeredgecolor','k')
set(hp(2),'markerfacecolor','y','markeredgecolor','k')
set(hp(3),'markerfacecolor','m','markeredgecolor','k')
set(hp(4),'markerfacecolor','c','markeredgecolor','k')
set(hp(5:7),'linewidth',2)
grid
title(['Equatorial Density ',titl],'fontsize',14,'fontweight','bold')
xlabel('L-shell','fontsize',12,'fontweight','bold')
ylabel('Log_{10} Electron density (m^{-3})','fontsize',12,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold')
legend(hp,'2 MLT','8 MLT','14 MLT','20 MLT','C. & A. saturated','C. & A. plasmatrough (0 MLT)   ','Gallagher 2000',1)

poinan=heden(:,edi1)>0;
ionrat1=heden(poinan,edi1)./hyden(poinan,edi1);
ionrat2=heden(poinan,edi2)./hyden(poinan,edi2);
ionrat3=heden(poinan,edi3)./hyden(poinan,edi3);
ionrat4=heden(poinan,edi4)./hyden(poinan,edi4);

if hes==1
 figure(12)
 hp=semilogy(lshm(poinan,1),ionrat1,'gs-',lshm(poinan,1),ionrat2,'yo-',lshm(poinan,1),ionrat3,'mp-',lshm(poinan,1),ionrat4,'c^-');
 set(hp(1),'markerfacecolor','g','markeredgecolor','k')
 set(hp(2),'markerfacecolor','y','markeredgecolor','k')
 set(hp(3),'markerfacecolor','m','markeredgecolor','k')
 set(hp(4),'markerfacecolor','c','markeredgecolor','k')
 grid
 title(['Equatorial He^+/H^+ Ratio ',titl],'fontsize',14,'fontweight','bold')
 xlabel('L-shell','fontsize',12,'fontweight','bold')
 ylabel('Log_{10} He^+/H^+ Ratio','fontsize',12,'fontweight','bold')
 set(gca,'fontsize',12,'fontweight','bold')
 legend(hp,'2 MLT','8 MLT','14 MLT','20 MLT',1)
end

figure(13)
surf(ones(size(dismat,1),1)*mltv,lshm,log10(eden));
title(['Equatorial Density ',titl],'fontsize',16,'fontweight','bold')
ylabel('L','fontsize',16,'fontweight','bold')
xlabel('MLT (hrs)','fontsize',16,'fontweight','bold')
zlabel('Log_{10} n(e^-) (m^-3)','fontsize',16,'fontweight','bold')
set(gca,'fontsize',14,'fontweight','bold')
view(142.5,30)

