%% plasirg.m is a script file that calls plasirg1.m to produce a regional
%% TEC map, as well as other outputs.
%%
%%
%% 17/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 - plasirg 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(' ')

lon1=input('First geographic longitude of region of interest (-180 to 180 degrees)? ');
lon2=input('Last geographic longitude of region of interest (-180 to 180 degrees)? ');
lonnop=input('Number of points to break geographic longitude into? ');
lat1=input('First geographic latitude of region of interest (-90 to 90 degrees)? ');
lat2=input('Last geographic latitude of region of interest (-90 to 90 degrees)? ');
latnop=input('Number of points to break geographic latitude into? ');
radnop=input('Number of points 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);
else 
 month=input('Month of interest? ');
 day=di;
end 
nop=input('Number of points to break each field line into during calculation? ');

%Determine the maximum allowed altitude.
re=6371;                                %Earth radi (km)
dno=daynum(day,month,year);
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.)? ');

lon1=deg2rad(lon1);
lon2=deg2rad(lon2);
lat1=deg2rad(lat1);
lat2=deg2rad(lat2);

[tec,tectot,mih,hytrh,hytrd,hetrh,hetrd,tecoh,techyp,techep,lonmat,latmat,eden,oxden,hyden,heden,disv,hes]=plasirg1(tecalt,lon1,lon2,lonnop,lat1,lat2,latnop,altmin,altmax,radnop,timmat,day,month,year,lshcp,nop);

timmat2=round(10*timmat)/10;
lonc=0.5*(lon1+lon2);
timl=timut2l(timmat,lonc);
timl=round(10*timl)/10;
lonc=rad2deg(lonc);

load coast
[lat,long]=maptriml(lat,long,rad2deg([lat1 lat2]),rad2deg([lon1 lon2]));

figure(1)
tecmin=floor(min(min(tectot)));
tecmax=ceil(max(max(tectot)));
tecv=tecmin: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
grid
xlabel('Geographic Longitude','fontsize',16,'fontweight','bold')
ylabel('Geographic Latitude','fontsize',16,'fontweight','bold')
titl=['UT=',num2str(timmat2),' (LT=',num2str(timl),' at Lon=',num2str(lonc),'^o) ',int2str(day),'/'int2str(month),'/',int2str(year)];
title(['TEC (TECU) ',titl],'fontsize',16,'fontweight','bold')
set(gca,'fontsize',14,'fontweight','bold','layer','top')

figure(2)
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];
[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
grid
xlabel('Geographic Longitude (degree)','fontsize',12,'fontweight','bold')
ylabel('Geographic Latitude (degree)','fontsize',12,'fontweight','bold')
title(['Mean Ionospheric Height (km) ',titl],'fontsize',15,'fontweight','bold')
set(gca,'fontsize',12,'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
grid
xlabel('Geographic Longitude (degree)','fontsize',12,'fontweight','bold')
ylabel('Geographic Latitude (degree)','fontsize',12,'fontweight','bold')
title(['O^+/H^+ Transition Height (km) ',titl],'fontsize',15,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold','layer','top')

if hes==1
 figure(4)
 hetrhmin=25*floor(min(min(hetrh/25)));
 hetrhmax=25*ceil(max(max(hetrh/25)));
 hetrhv=[100:25:450 500:50:800 900 1000:200:1800 2000:500:5000];
 hetrhv(hetrhv<hetrhmin | hetrhv>hetrhmax)=[];
 hetrhv=[hetrhmin hetrhv];
 [c,h]=contourf(rad2deg(lonmat),rad2deg(latmat),hetrh,hetrhv);
 hold on
 h1=plot(long,lat,'w');
 set(h1,'linewidth',1.5)
 hold off
 clabel(c,h)
 colorbar
 grid
 xlabel('Geographic Longitude (degree)','fontsize',12,'fontweight','bold')
 ylabel('Geographic Latitude (degree)','fontsize',12,'fontweight','bold')
 title(['O^+/He^+ Transition Height (km) ',titl],'fontsize',15,'fontweight','bold')
 set(gca,'fontsize',12,'fontweight','bold','layer','top')
end

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
grid
xlabel('Geographic Longitude (degree)','fontsize',12,'fontweight','bold')
ylabel('Geographic Latitude (degree)','fontsize',12,'fontweight','bold')
title(['TEC above O^+/H^+ Transition Height (TECU) ',titl],'fontsize',15,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold','layer','top')

figure(6)
tecmax=ceil(max(max(techyp)));
tecv2=0:2.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
grid
xlabel('Geographic Longitude (degree)','fontsize',12,'fontweight','bold')
ylabel('Geographic Latitude (degree)','fontsize',12,'fontweight','bold')
title(['Percentage of TEC due to H^+ ',titl],'fontsize',15,'fontweight','bold')
set(gca,'fontsize',12,'fontweight','bold','layer','top')

if hes==1
 figure(7)
 tecmax=ceil(max(max(techep)));
 tecv2=0:0.2: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
 grid
 xlabel('Geographic Longitude (degree)','fontsize',12,'fontweight','bold')
 ylabel('Geographic Latitude (degree)','fontsize',12,'fontweight','bold')
 title(['Percentage of TEC due to He^+ ',titl],'fontsize',15,'fontweight','bold')
 set(gca,'fontsize',12,'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
  grid
  xlabel('Geographic Longitude (degree)','fontsize',12,'fontweight','bold')
  ylabel('Geographic Latitude (degree)','fontsize',12,'fontweight','bold')
  title(['TEC (TECU) - ',num2str(tecalt2(i)),' km to ',num2str(tecalt2(i+1)),' km'],'fontsize',15,'fontweight','bold')
  set(gca,'fontsize',12,'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 

