%% plasisc.m is a script file that calls plasisc1.m to produce a H+, He+, O+
%% and e- density magnetic longitude slice, as well as other plots.
%% 
%% 16/11/2001 - updated 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 all
close all

lonlsh=input('Magnetic longitude of interest (degrees)? ');
latmnop=input('Number of latitudinal divisions to break -90 to 90 into? ');
radnop=input('Number of divisions to break radial ray paths into? ');
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 
nop=input('Number of points to break each field line into during calculation? ');

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

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

lonlsh=deg2rad(lonlsh);

[eden,oxden,hyden,heden,latmatm,dismat,hes]=plasisc1(lonlsh,latmnop,altmin,altmax,radnop,timmat,day,month,year,nop);

kp=solarval(timmat,dno,year);
kp=round(10*kp)/10;
[xmat,ymat]=pol2cart(latmatm,dismat);
cf=7.5/18.5;
re=6371;
lsx=max(max(dismat))/re;
tecalt=tecalt(:)';
[tec,mih,tectot]=teccal1(eden,dismat(1,:),tecalt);
hytec=teccal1(hyden,dismat(1,:),0);
[hytrh,hytrd]=tranht(dismat-re,oxden,hyden);
if hes==1
 hetec=teccal1(heden,dismat(1,:),0);
 [hetrh,hetrd]=tranht(dismat-re,oxden,heden);
else
 hetec=nan*ones(size(hytec));
 hetrh=nan*ones(size(hytrh));
 hetrd=nan*ones(size(hytrd));
end

tectr=teccal2(eden,dismat(1,:),hytrh);
timm=timut2ml2(timmat,lonlsh,dno);
timm=round(timm*100)/100;
timmat2=round(10*timmat)/10;

svm=ones(size(dismat));
[lonmat,latmat]=mag2geo3(lonlsh*svm,latmatm);
[tmpi,tmpe]=temppaw(lonmat,latmat,dismat,timmat*svm,year*svm,month*svm,day*svm);
tmpe((dismat-re)<350)=nan;

figure(1)
pcolor(xmat/re,ymat/re,log10(eden))
shading flat
axis([0 lsx -cf*lsx cf*lsx]) 
colorbar
hold on
[xme,yme]=pol2cart(latmatm(:,1),1);
fill(xme,yme,'g')
hold off
xlabel('Earth radii (L-shell)','fontsize',12,'fontweight','bold')
ylabel('Earth radii (L-shell)','fontsize',12,'fontweight','bold')
titl=['UT=',num2str(timmat2),' (MLT=',num2str(timm),')  ',int2str(day),'/',int2str(month),'/',int2str(year),' Lonm=',num2str(rad2deg(lonlsh)),'^{o} Kp=',num2str(kp)];
title(['Log_{10} e^- den.  ',titl],'fontsize',12,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold')

latgps=deg2rad(-55:5:55);
[xvgps,yvgps]=pol2cart(latgps,(re+20200)*ones(size(latgps)));

figure(2)
pcolor(xmat/re,ymat/re,log10(eden))
shading flat
hold on
h21=plot([0 max(max(xmat/re))],[0 0],'color','k','linewidth',2);
h22=plot(xvgps/re,yvgps/re,'ok','markerfacecolor','w');
[xme,yme]=pol2cart(latmatm(:,1),1);
fill(xme,yme,'g')
text(0.2,0.3,'Earth','fontsize',14,'color','k','fontweight','bold')
text(1,-0.5,'\leftarrow Ionosphere','fontsize',14,'color','k','fontweight','bold')
text(2,0.5,'Plasmasphere','fontsize',14,'color','k','fontweight','bold')
hold off
axis([0 lsx -cf*lsx cf*lsx]) 
colorbar
xlabel('Earth radii','fontsize',16,'fontweight','bold')
ylabel('Earth radii','fontsize',16,'fontweight','bold')
titl=['UT=',num2str(timmat2),' (MLT=',num2str(timm),')  ',int2str(day),'/',int2str(month),'/',int2str(year),' Lonm=',num2str(rad2deg(lonlsh)),'^{o} Kp=',num2str(kp)];
title(['Log_{10} e^- den.  ',titl],'fontsize',16,'fontweight','bold')
set(gca,'fontsize',14,'fontweight','bold')
legend([h21 h22],'Mag. Equator','GPS orbital plane',4)

oxden2=oxden;
oxden2(oxden2==0)=NaN;
figure(3)
pcolor(xmat/re,ymat/re,log10(oxden2))
shading flat
axis([0 lsx -cf*lsx cf*lsx]) 
colorbar
hold on
[xme,yme]=pol2cart(latmatm(:,1),1);
fill(xme,yme,'g')
hold off
xlabel('Earth radii (L-shell)','fontsize',12,'fontweight','bold')
ylabel('Earth radii (L-shell)','fontsize',12,'fontweight','bold')
title(['Log_{10} O^+ den.  ',titl],'fontsize',12,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold')

hyden2=hyden;
hyden2(hyden2==0)=NaN;
figure(4)
pcolor(xmat/re,ymat/re,log10(hyden2))
shading flat
axis([0 lsx -cf*lsx cf*lsx]) 
colorbar
hold on
[xme,yme]=pol2cart(latmatm(:,1),1);
fill(xme,yme,'g')
hold off
xlabel('Earth radii (L-shell)','fontsize',12,'fontweight','bold')
ylabel('Earth radii (L-shell)','fontsize',12,'fontweight','bold')
title(['Log_{10} H^+ den.  ',titl],'fontsize',12,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold')

if hes==1
 heden2=heden;
 heden2(heden2==0)=NaN;
 figure(5)
 pcolor(xmat/re,ymat/re,log10(heden2))
 shading flat
 axis([0 lsx -cf*lsx cf*lsx]) 
 colorbar
 hold on
 [xme,yme]=pol2cart(latmatm(:,1),1);
 fill(xme,yme,'g')
 hold off
 xlabel('Earth radii (L-shell)','fontsize',12,'fontweight','bold')
 ylabel('Earth radii (L-shell)','fontsize',12,'fontweight','bold')
 title(['Log_{10} He^+ den.  ',titl],'fontsize',12,'fontweight','bold')
 set(gca,'fontsize',12,'fontweight','bold')
end

cv=['r';'g';'b';'c';'m';'y'];
cv2=cv(1:(length(tecalt)+1));
cv2=flipud(cv2);
latd=rad2deg(latmatm(:,1));
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(6)
h0=area(latd,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([-90 90 0 10*ceil(max(tectot)/10)])
xlabel('Magnetic Latitude','fontsize',12,'fontweight','bold')
ylabel('Total Electron Content (TECU)','fontsize',12,'fontweight','bold')
title(['TEC Altitude bands: ',titl],'fontsize',12,'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(7)
h0=area(latd,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([-90 90 0 10*ceil(max(tectot)/10)])
xlabel('Magnetic Latitude','fontsize',16,'fontweight','bold')
ylabel('Total Electron Content (TECU)','fontsize',16,'fontweight','bold')
title(['TEC Altitude bands: ',titl],'fontsize',16,'fontweight','bold')
set(gca,'fontsize',14,'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(8)
hp=plot(latd,mih,'bo-');
grid
xlabel('Magnetic Latitude','fontsize',12,'fontweight','bold')
ylabel('Median ionospheric height (km)','fontsize',12,'fontweight','bold')
title(['Median Ionospheric Height: ',titl],'fontsize',12,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold')
set(hp,'linewidth',2)
set(hp(1),'markerfacecolor','b')

figure(9)
subplot(1,2,1)
hp=plot(latd,hytrh,'bo-',latd,hetrh,'^r-');
grid
title('Transition Height','fontsize',12,'fontweight','bold')
xlabel('Magnetic Latitude','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')
legend(hp,'O^+/H^+','O^+/He^+',1)
subplot(1,2,2)
hytrd(hytrd==0)=NaN;
hp=plot(latd,log10(hytrd),'bo-',latd,log10(hetrd),'r^-');
grid
xlabel('Magnetic Latitude','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')
legend(hp,'O^+/H^+','O^+/He^+',1)

figure(10)
subplot(2,1,1)
h0=area(latd,tectr);
set(h0(2),'facecolor','c')
set(h0(1),'facecolor','r')
grid
axis([-90 90 0 10*ceil(max(tectot)/10)])
xlabel('Magnetic Latitude','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',1)
subplot(2,1,2)
h0=plot(latd,tectr(:,2),'b^-');
grid
axis([-90 90 0 10*ceil(max(tectr(:,2))/10)])
xlabel('Magnetic Latitude','fontsize',12,'fontweight','bold')
ylabel('Total Electron Content (TECU)','fontsize',12,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold')
set(h0,'linewidth',2,'markerfacecolor','b')
legend(h0,'TEC above',1)

figure(11)
h0=plot(latd,100*hytec./tectot,'r^-',latd,100*hetec./tectot,'bo-');
grid
axis([-90 90 0 10*ceil(100*max(max(hytec./tectot),max(hetec./tectot))/10)])
title('Percentage of TEC due to','fontsize',12,'fontweight','bold')
xlabel('Magnetic Latitude','fontsize',12,'fontweight','bold')
ylabel('% of TEC','fontsize',12,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold')
set(h0(1),'linewidth',2,'markerfacecolor','r')
set(h0(2),'linewidth',2,'markerfacecolor','b')
legend(h0,'H^+','He^+',1)

poseq=round(latmnop/2);       
oxeq=oxden2(poseq,:);
hyeq=hyden2(poseq,:);
if hes==1
 heeq=heden2(poseq,:);
else
 heeq=nan*ones(size(hyeq));
end
oxpol=oxden2(1,:);
hypol=hyden2(1,:);
if hes==1
 hepol=heden2(1,:);
else
 hepol=nan*ones(size(hypol));
end
alt=dismat(1,:)-re;
figure(12)
subplot(1,2,1)
hp=semilogx(oxeq,alt,'bo-',hyeq,alt,'r^-',heeq,alt,'gs-');
grid
title('Magnetic Equator','fontsize',12,'fontweight','bold')
xlabel('Ion Density (m^{-3})','fontsize',10,'fontweight','bold')
ylabel('Altitude (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')
set(hp(3),'markerfacecolor','g')
legend(hp,'O^+','H^+','He^+')
subplot(1,2,2)
hp=semilogx(oxpol,alt,'bo-',hypol,alt,'r^-',hepol,alt,'gs-');
grid
title('Northern Magnetic Pole','fontsize',12,'fontweight','bold')
xlabel('Ion Density (m^{-3})','fontsize',10,'fontweight','bold')
ylabel('Altitude (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')
set(hp(3),'markerfacecolor','g')
legend(hp,'O^+','H^+','He^+')

svy=ones(1,size(dismat,2));
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;

figure(13)
hp=semilogy(lshm(1,:),hyeq,'b',lshm(1,:),heeq,'g',lshm(1,:),(oxeq+hyeq+heeq),'r',lshm(1,:),1e6*eca,'sb',lshm(1,:),1e6*eptca,'^m',lshm(1,:),1e6*edeng,'og');
hpl=length(hp);
set(hp(1:3),'linewidth',2)
set(hp(4),'markerfacecolor','b')
set(hp(5),'markerfacecolor','m')
set(hp(6),'markerfacecolor','g')
grid
title(['Equatorial Density ',titl],'fontsize',12,'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,'GPID H^+','GPID He^+','GPID e^- profiles','C. & A. saturated','C. & A. plasmatrough (0 MLT)   ','Gallagher 1999',1)

figure(14)
pcolor(xmat/re,ymat/re,tmpe)
shading flat
axis([0 lsx -cf*lsx cf*lsx]) 
colorbar
hold on
[xme,yme]=pol2cart(latmatm(:,1),1);
fill(xme,yme,'g')
hold off
xlabel('Earth radii (L-shell)','fontsize',12,'fontweight','bold')
ylabel('Earth radii (L-shell)','fontsize',12,'fontweight','bold')
titl=['UT=',num2str(timmat2),' (MLT=',num2str(timm),')  ',int2str(day),'/',int2str(month),'/',int2str(year),' ',num2str(rad2deg(lonlsh)),'^{o} mag. long Kp=',num2str(kp)];
title(['e^- temperature (K)  ',titl],'fontsize',12,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold')


