%% plasivert.m calculates the required densities for a vertical
%% ray path above a given geographic location for the given day.

clear all
close all

tims=input('Time step to be used (fractional hours)? ');
long=input('Geographic longitude (degrees)? ');
latg=input('Geographic latitude (degrees)? ');
noprp=input('Number of points to break ray paths into? ');
dno=input('Day number? ');
year=input('Year? ');
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.)? ');


timmat=0:tims:(24-tims);
sv=ones(size(timmat));
loni=deg2rad(long)*sv;
lati=deg2rad(latg)*sv;
alti=altmin*sv;
altf=altmax*sv;
lt=0;
[day,month]=daynum2(dno,year);

[tec,mih,lonmih,latmih,eden,oxden,hyden,heden,lonmat,latmat,dismat]=plasiray1(loni,loni,lati,lati,alti,altf,noprp,lt,timmat,day*sv,month*sv,year*sv,nop);

re=6371;
[tec,mih,tectot]=teccal1(eden,dismat(1,:),tecalt);
[hytrh,hytrd]=tranht(dismat-re,oxden,hyden);
[hetrh,hetrd]=tranht(dismat-re,oxden,heden);

timmat=[timmat 24];
tec=[tec;tec(1,:)];
tectot=[tectot;tectot(1,:)];
mih=[mih;mih(1)];
hytrh=[hytrh;hytrh(1)];
hytrd=[hytrd;hytrd(1)];
hetrh=[hetrh;hetrh(1)];
hetrd=[hetrd;hetrd(1)];

timl=timut2l(0,loni(1));
if timl<=12
 timl=round(10*timl)/10;
 xlab=['UT (hours) (LT = UT + ',num2str(timl),')'];
else 
 timl=timl-24;
 timl=abs(round(10*timl)/10);
 xlab=['UT (hours) (LT = UT - ',num2str(timl),')'];
end 
titl=[int2str(day),'/'int2str(month),'/',int2str(year),'   Lon=',num2str(long),'^o Lat=',num2str(latg),'^o'];

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(1)
h0=area(timmat,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(xlab,'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')
set(gca,'xtick',0:4:24)
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(2)
h0=area(timmat,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(xlab,'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')
set(gca,'xtick',0:4:24)
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(3)
hp=plot(timmat,mih,'bo-');
grid
xlabel(xlab,'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')
set(gca,'xtick',0:4:24)

figure(4)
subplot(1,2,1)
hp=plot(timmat,hytrh,'bo-',timmat,hetrh,'r^-');
grid
xlabel(xlab,'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(timmat,log10(hytrd),'bo-',timmat,log10(hetrd),'^r-');
grid
xlabel(xlab,'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')
set(gca,'xtick',0:4:24)
legend(hp,'O^+/H^+','O^+/He^+',1) 

