%Demonstrates the difference between vertical and slant TEC along a latitude slice
%for two "observation stations" at 30 and 24 deg. 

clear all 
close all

lats=-deg2rad(15:40);
latsd=rad2deg(lats);
sv=ones(size(lats));
lons=deg2rad(-70)*sv;

timmat=7;
dno=228;
year=1995;
[day,month]=daynum2(dno,year);

noprp=50;
lt=0;
iricf=1;
nop=50;

loni=lons;
lonf=lons;
lati=lats;
latf=lats;
alti=0*sv;
altf=20200*sv;

re=6371;

altsh=400;

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

%Calculate slant TEC where the ray paths from the observing stations are assumed to pass through 400km 
%when the vertical TEC is calculated.
[rg30,azs30,els30]=satloc2(timmat*sv,day*sv,month*sv,year*sv,lons,lats,(re+altsh)*sv,lons,deg2rad(-30)*sv,re*sv);
[lons30,lats30]=satloc(rad2deg(lons),(-30)*sv,rad2deg(azs30),rad2deg(els30),20200*sv);
[tec30,mih30]=plasiray1(lons,deg2rad(lons30),deg2rad(-30)*sv,deg2rad(lats30),0*sv,20200*sv,noprp,lt,timmat*sv,day*sv,month*sv,year*sv,iricf,nop);
sf30=1./sqrt(1-(cos(els30)./(1+(altsh/re))).^2);
tec30v=(tec30')./sf30;

[rg24,azs24,els24]=satloc2(timmat*sv,day*sv,month*sv,year*sv,lons,lats,(re+altsh)*sv,lons,deg2rad(-24)*sv,re*sv);
[lons24,lats24]=satloc(rad2deg(lons),(-24)*sv,rad2deg(azs24),rad2deg(els24),20200*sv);
[tec24,mih24]=plasiray1(lons,deg2rad(lons24),deg2rad(-24)*sv,deg2rad(lats24),0*sv,20200*sv,noprp,lt,timmat*sv,day*sv,month*sv,year*sv,iricf,nop);
sf24=1./sqrt(1-(cos(els24)./(1+(altsh/re))).^2);
tec24v=(tec24')./sf24;

timlt=timut2l(timmat,lons(1));

figure(1)
subplot(2,1,1)
h1=plot(latsd,mih,'o-b');
set(h1,'linewidth',2,'markerfacecolor','b')
grid
axis([min(latsd) max(latsd) 50*floor(min(mih)/50) 50*ceil(max(mih)/50)])
set(gca,'fontsize',12,'fontweight','bold')
xlabel('Geographic latitude (deg)','fontsize',12,'fontweight','bold')
ylabel('Mean Ionospheric Height (km)','fontsize',12,'fontweight','bold')
tits=['Lon=',int2str(rad2deg(lons(1))),'^o   LT=',num2str(timlt,2),'   UT=',num2str(timmat,2),'   ',int2str(day),'/',int2str(month),'/',int2str(year)];
title(tits,'fontsize',14,'fontweight','bold')

subplot(2,1,2)
h2=plot(latsd,tec,'o-b',latsd,tec30v,'^-r',latsd,tec24v,'d-g',-30,0,'sr',-24,0,'sg');
set(h2(1),'linewidth',2,'markerfacecolor','b')
set(h2(2),'linewidth',2,'markerfacecolor','r')
set(h2(3),'linewidth',2,'markerfacecolor','g')
set(h2(4),'markersize',10,'markerfacecolor','r')
set(h2(5),'markersize',10,'markerfacecolor','g')
grid
axis([min(latsd) max(latsd) 0 2*ceil(max([tec' tec30v tec24v])/2)])
set(gca,'fontsize',12,'fontweight','bold')
xlabel('Geographic latitude (deg)','fontsize',12,'fontweight','bold')
ylabel('TEC (TECU)','fontsize',12,'fontweight','bold')
title(['Ionospheric Shell Height ',int2str(altsh),' km'],'fontsize',14,'fontweight','bold')
legend('True vertical','Slant to vertical 30^o station','Slant to vertical 24^o station',4)

