%% plasigb.m is a script file that calls plasigb1.m to produce a global
%% TEC map and other outputs.
%%
%% 20/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

disp(' ')
disp(' ')
disp('NOTE - plasigb uses a standard size longitude and latitude steps,')
disp('the values inputted by the user are used to interpolate the densities')
disp('from the standard grid.  Consequently computation time basically ')
disp('independent of grid size given by user.')
disp(' ')
disp(' ')

lonss=input('Geographic longitude step size (degrees - must divide exactly into 360)? ');
latss=input('Geographic latitude step size (degrees - must divide exactly into 180)? ');
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.)? ');

[tec,tectot,mih,hytrh,hytrd,hetrh,hetrd,tecoh,techyp,techep,lonmat,latmat,eden,oxden,hyden,heden,disv,hes]=plasigb1(tecalt,lonss,latss,altmin,altmax,radnop,timmat,day,month,year,lshcp,nop);

timmat2=round(10*timmat)/10;

load coast

figure(1)
tecmin=floor(min(min(tectot)));
tecmax=ceil(max(max(tectot)));
tecv=tecmin:2:tecmax;
[c,h]=contourf(rad2deg(lonmat),rad2deg(latmat),tectot,tecv);
hold on
h1=plot(long,lat,'w');
set(h1,'linewidth',1.5)
hold off
clabel(c,h)
colorbar('horiz')
grid
xlabel('Geographic Longitude','fontsize',14,'fontweight','bold')
ylabel('Geographic Latitude','fontsize',14,'fontweight','bold')
titl=['UT=',num2str(timmat2),' ',int2str(day),'/',int2str(month),'/',int2str(year)];
title(['Total Vertical TEC (TECU) ',titl],'fontsize',16,'fontweight','bold')
set(gca,'fontsize',14,'fontweight','bold','layer','top')

mihmin=25*floor(min(min(mih/25)));
mihmax=25*ceil(max(max(mih/25)));
mihv=[100:25:450 500:50:800 900 1000:200:1800 2000:500:5000];
mihv(mihv<mihmin | mihv>mihmax)=[];
mihv=[mihmin mihv];
figure(2)
[c,h]=contourf(rad2deg(lonmat),rad2deg(latmat),mih,mihv);
hold on
h1=plot(long,lat,'w');
set(h1,'linewidth',1.5)
hold off
clabel(c,h)
colorbar('horiz')
grid
xlabel('Geographic Longitude','fontsize',14,'fontweight','bold')
ylabel('Geographic Latitude','fontsize',14,'fontweight','bold')
title(['Mean Ionospheric Height (km) ',titl],'fontsize',16,'fontweight','bold')
set(gca,'fontsize',14,'fontweight','bold','layer','top')

figure(3)
hytrhmin=25*floor(min(min(hytrh/25)));
hytrhmax=25*ceil(max(max(hytrh/25)));
hytrhv=[100:25:450 500:50:800 900 1000:200:1800 2000:500:5000];
hytrhv(hytrhv<hytrhmin | hytrhv>hytrhmax)=[];
hytrhv=[hytrhmin hytrhv];
[c,h]=contourf(rad2deg(lonmat),rad2deg(latmat),hytrh,hytrhv);
hold on
h1=plot(long,lat,'w');
set(h1,'linewidth',1.5)
hold off
clabel(c,h)
colorbar('horiz')
grid
xlabel('Geographic Longitude','fontsize',14,'fontweight','bold')
ylabel('Geographic Latitude','fontsize',14,'fontweight','bold')
title(['O^+/H^+ Transition Height (km) ',titl],'fontsize',16,'fontweight','bold')
set(gca,'fontsize',14,'fontweight','bold','layer','top')

figure(4)
[c,h]=contourf(rad2deg(lonmat),rad2deg(latmat),log10(hytrd),15);
hold on
h1=plot(long,lat,'w');
set(h1,'linewidth',1.5)
hold off
clabel(c,h)
colorbar('horiz')
grid
xlabel('Geographic Longitude','fontsize',14,'fontweight','bold')
ylabel('Geographic Latitude','fontsize',14,'fontweight','bold')
title(['Log_{10} O^+/H^+ Transition Density (m^{-3}) ',titl],'fontsize',16,'fontweight','bold')
set(gca,'fontsize',14,'fontweight','bold','layer','top')

figure(5)
tecoh2=squeeze(tecoh(:,:,2));
tecmin=floor(min(min(tecoh2)));
tecmax=ceil(max(max(tecoh2)));
tecv2=tecmin:0.5:tecmax;
[c,h]=contourf(rad2deg(lonmat),rad2deg(latmat),tecoh2,tecv2);
hold on
h1=plot(long,lat,'w');
set(h1,'linewidth',1.5)
hold off
clabel(c,h)
colorbar('horiz')
grid
xlabel('Geographic Longitude','fontsize',14,'fontweight','bold')
ylabel('Geographic Latitude','fontsize',14,'fontweight','bold')
title(['TEC above O^+/H^+ Transition Height (TECU) ',titl],'fontsize',16,'fontweight','bold')
set(gca,'fontsize',14,'fontweight','bold','layer','top')

figure(6)
tecmax=ceil(max(max(techyp)));
tecv2=0:5:tecmax;
[c,h]=contourf(rad2deg(lonmat),rad2deg(latmat),techyp,tecv2);
hold on
h1=plot(long,lat,'w');
set(h1,'linewidth',1.5)
hold off
clabel(c,h)
colorbar('horiz')
grid
xlabel('Geographic Longitude','fontsize',14,'fontweight','bold')
ylabel('Geographic Latitude','fontsize',14,'fontweight','bold')
title(['Percentage of TEC due to H^+ ',titl],'fontsize',16,'fontweight','bold')
set(gca,'fontsize',14,'fontweight','bold','layer','top')

if hes==1
 figure(7)
 tecmax=ceil(max(max(techep)));
 tecv2=0:1:tecmax;
 [c,h]=contourf(rad2deg(lonmat),rad2deg(latmat),techep,tecv2);
 hold on
 h1=plot(long,lat,'w');
 set(h1,'linewidth',1.5)
 hold off
 clabel(c,h)
 colorbar('horiz')
 grid
 xlabel('Geographic Longitude','fontsize',14,'fontweight','bold')
 ylabel('Geographic Latitude','fontsize',14,'fontweight','bold')
 title(['Percentage of TEC due to He^+ ',titl],'fontsize',16,'fontweight','bold')
 set(gca,'fontsize',14,'fontweight','bold','layer','top')
end

if tecalt~=0
 tecalt2=[0 tecalt altmax];
 for i=1:(length(tecalt)+1)
  figure(i+7)
  teci=squeeze(tec(:,:,i));
  [c,h]=contourf(rad2deg(lonmat),rad2deg(latmat),teci);
  hold on
  h1=plot(long,lat,'w');
  set(h1,'linewidth',1.5)
  hold off
  clabel(c,h)
  colorbar('horiz')
  grid
  xlabel('Geographic Longitude','fontsize',14,'fontweight','bold')
  ylabel('Geographic Latitude','fontsize',14,'fontweight','bold')
  title([titl,'  TEC (TECU) - ',num2str(tecalt2(i)),' km to ',num2str(tecalt2(i+1)),' km'],'fontsize',16,'fontweight','bold')
  set(gca,'fontsize',14,'fontweight','bold','layer','top')
 end 
end 

[lonmatm,latmatm]=geo2mag3(lonmat,latmat);

if tecalt~=0
 tecalt2=[0 tecalt altmax];
 for k=1:(length(tecalt)+1)
  figure(i+7+k)
  teci=squeeze(tec(:,:,k));
  plot(rad2deg(latmatm),teci,'ro');
  grid
  xlabel('Magnetic Latgitude (degree)','fontsize',12,'fontweight','bold')
  ylabel('TEC (TECU)','fontsize',12,'fontweight','bold')
  title([titl,' Global TEC - ',num2str(tecalt2(i)),' km to ',num2str(tecalt2(i+1)),' km'],'fontsize',15,'fontweight','bold')
  set(gca,'fontsize',12,'fontweight','bold','layer','top')
 end 
end 
