%% sacccp1.m provides a comparison between SAC-C TEC data and
%% GPID prediction. If the GPID prediction has been previously
%% calculated, it is loaded in, otherwise plasiray.m is used
%% to calculate the required values.
%%

close all
clear all

%Date vectors
%dnov=[190];
%yearv=[2001];
%dnoi=menu('Day number of interest','190');
%dno=dnov(dnoi);
%yeari=menu('Year of interest','2002');
%year=yearv(yeari);

dno=190;
year=2001;

dats=menu('Use which data set','Full','Reduced');
dats=dats-1;

%Load in data from JPL file
[tec,lonoe,latoe,disoe,longp,latgp,disgp,timmat,day,month,year,svn,tecbc]=saccrd(dno,year,dats);

re=6371;                           %Earth radi (km)

%Determine if the required GPID results have already been calculated. If they
%have, load them in. Otherwise use plasiray to determine their value.
if dats==0
 strn=['tecres',int2str(dno),int2str(year(1))];
else
 strn=['tecres',int2str(dno),int2str(year(1)),'r'];
end
fe=feval('exist',['JPL\',strn,'.mat']);

if fe==2
 eval(['load JPL\',strn])
else
 dp=length(tec);
 noprp=50;                    %number of points to break each radial path
 lt=0;                        %use log spacing
 nop=50;                      %number of points for each field line
 nopl=40;                     %number of ray to simulate per call to plasiray
 nol=ceil(dp/nopl);           %number of calls to plasiray required
 
 %Result vectors
 tecgp=zeros(dp,1);
 mih=tecgp;
 lonmih=tecgp;
 latmih=tecgp;
 eden=zeros(dp,noprp);
 oxden=eden;
 hyden=eden;
 heden=eden;
 lonmat=eden;
 latmat=eden;
 dismat=eden;
 
 tic
 for j=1:nol
  [j nol]
  poil=(j-1)*nopl+1;          %lower index required
  if j~=nol
   poiu=poil+nopl-1;          %upper index required
  else
   poiu=dp;
  end
  poi=poil:poiu;              %required index for j'th loop
  loni=lonoe(poi);
  lonf=longp(poi);
  lati=latoe(poi);
  latf=latgp(poi);
  alti=disoe(poi)-re;
  altf=disgp(poi)-re;
  timmat2=timmat(poi);
  day2=day(poi);
  month2=month(poi);
  year2=year(poi);
  
  %Calculate required parameters
  [tecgpj,mihj,lonmihj,latmihj,edenj,oxdenj,hydenj,hedenj,lonmatj,latmatj,dismatj]=plasiray1(loni,lonf,lati,latf,alti,altf,noprp,lt,timmat2,day2,month2,year2,nop);
 
  tecgp(poi)=tecgpj;
  mih(poi)=mihj;
  lonmih(poi)=lonmihj;
  latmih(poi)=latmihj;
  eden(poi,:)=edenj;
  oxden(poi,:)=oxdenj;
  hyden(poi,:)=hydenj;
  heden(poi,:)=hedenj;
  lonmat(poi,:)=lonmatj;
  latmat(poi,:)=latmatj;
  dismat(poi,:)=dismatj;
  
  toc
  
 end
 
 eval(['save JPL\',strn,' dno nop tecgp mih lonmih latmih eden oxden hyden heden lonmat latmat dismat'])
 
end

%Calculate the range, azimuth and elevation of the GPS satellite relative to SAC-C.
[rg,azs,els]=satloc2(timmat,day,month,year,lonmat(:,nop),latmat(:,nop),dismat(:,nop),lonmat(:,1),latmat(:,1),dismat(:,1));

tecgp(tecgp==0)=nan;               %replace GPID zero values with NaN

%Magnetic co-ords
[lonm,latm]=geo2mag3(lonmat(:,1),latmat(:,1));
latmd=rad2deg(latm);
timl=timut2l(timmat,lonmat(:,1));
elsd=rad2deg(els);

ax1=5*max(ceil(tecgp/5));
ax2=5*max(ceil(tec/5));
ax3=max(ax1,ax2);

men1='Total SAC-C TEC vs Total GPID TEC';
men2='Particular Elevation Mask';
men3='Elevations masks in 10 degree steps';
men4='Magnetic latitude variation';
men5='TEC diff vs LT';
 
nopi=menu('Plotting Option',men1,men2,men3,men4,men5);

if nopi==1

 %Conduct line fit to SAC-C vs GPID TECU
 [c1,c0,c1e,c0e]=linfit(tecgp,tec);
 tecgpf=linspace(0,ax1,100);
 tecf=polyval([c1 c0],tecgpf);
 
 %Save GPID TEC to ascii file
 %fid=fopen('paw_tec.txt','w');
 %fprintf(fid,'%6.2f\n',tecgp);
 %fclose(fid)
 
 poid=timl>=6 & timl<18;
 poin=timl>=18 | timl<6;
  
 figure(1)
 h1=plot(tecgp,tec,'.b',tecgpf,tecf,'k','linewidth',2);
 grid
 axis([0 ax1 0 ax2])
 xlabel('GPID TEC (TECU)','fontweight','bold','fontsize',16)
 ylabel('SAC-C TEC (TECU)','fontweight','bold','fontsize',16)
 %title(['SAC-C vs GPID TEC (all data) (',iris,') Day=',int2str(dno),' 1999'],'fontweight','bold','fontsize',16)
 set(gca,'fontweight','bold','fontsize',14)
 fls=['gradient=',num2str(c1,'%6.3f'),'\pm',num2str(c1e,'%6.3f'),'   intercept=',num2str(c0,'%6.2f'),'\pm',num2str(c0e,'%6.2f')];
 legend(h1(2),fls,4)
 
 figure(2)
 subplot(2,2,1)
 h12=plot(tec(poid),latmd(poid),'.b',tec(poin),latmd(poin),'.r');
 grid 
 set(gca,'fontweight','bold')
 set(gca,'ytick',-90:30:90)
 xlabel('SAC-C TEC (TECU)','fontweight','bold')
 ylabel('Mag. Lat.','fontweight','bold')
 title(['SAC-C Dno=',int2str(dno),' 1999'],'fontweight','bold')
 axis([0 ax3 -90 90])
 legend(h12,'Day','Night',4)

 subplot(2,2,2)
 h12=plot(tecgp(poid),latmd(poid),'.b',tecgp(poin),latmd(poin),'.r');
 grid 
 set(gca,'fontweight','bold')
 set(gca,'ytick',-90:30:90)
 xlabel('GPID TEC (TECU)','fontweight','bold')
 ylabel('Mag. Lat.','fontweight','bold')
 title(['GPID (',iris,') Dno=',int2str(dno),' 1999'],'fontweight','bold')
 axis([0 ax3 -90 90])
 legend(h12,'Day','Night',4)

 subplot(2,2,3)
 plot(tecgp,timmat,'.')
 grid
 set(gca,'ytick',0:4:24)
 set(gca,'fontweight','bold')
 xlabel('GPID TEC (TECU)','fontweight','bold')
 ylabel('Time (UT)','fontweight','bold')
 title(['UT (',iris,') Dno=',int2str(dno),' 1999'],'fontweight','bold')
 axis([0 ax1 0 24])

 subplot(2,2,4)
 plot(tecgp,timl,'.')
 grid
 set(gca,'ytick',0:4:24)
 set(gca,'fontweight','bold')
 xlabel('GPID TEC (TECU)','fontweight','bold')
 ylabel('Time (LT)','fontweight','bold')
 title(['LT (SAC-C) (',iris,') Dno=',int2str(dno),' 1999'],'fontweight','bold')
 axis([0 ax1 0 24])


elseif nopi==2
 
 dgv=input('Minimum elevation angle to use from SAC-C (degrees)? ');
 
 poime=elsd>=dgv;
 poid=poime & (timl>=9 & timl<=15);
 poin=poime & (timl>=21 | timl<=3);
 
 ax1d=5*max(ceil(tecgp(poid)/5));
 ax1n=5*max(ceil(tecgp(poid)/5));
 ax2d=5*max(ceil(tec(poid)/5));
 ax2n=5*max(ceil(tec(poin)/5));
 ax1x=max(ax1d,ax1n);
 ax2x=max(ax2d,ax2n);
 ax3x=max(ax1x,ax2x);

 tecgpf=linspace(0,ax1,100);
  
 [c1d,c0d,c1ed,c0ed]=linfit(tecgp(poid),tec(poid));
 tecfd=polyval([c1d c0d],tecgpf);
 
 [c1n,c0n,c1en,c0en]=linfit(tecgp(poin),tec(poin));
 tecfn=polyval([c1n c0n],tecgpf);
 
 figure(1)
 h1=plot(tecgp(poid),tec(poid),'xb',tecgp(poin),tec(poin),'or',tecgpf,tecfd,'b',tecgpf,tecfn,'r--');
 grid
 set(gca,'xtick',0:2:ax1)
 set(gca,'ytick',0:2:ax2)
 set(gca,'fontweight','bold','fontsize',10)
 axis([0 ax1x 0 ax2x])
 xlabel('GPID TEC (TECU)','fontweight','bold','fontsize',12)
 ylabel('SAC-C TEC (TECU)','fontweight','bold','fontsize',12)
 title(['+',num2str(dgv,2),'^o SAC-C elevation mask angle (',iris,') Dno=',int2str(dno),' 1999'],'fontweight','bold','fontsize',14)
 days=['Day fit    grad=',num2str(c1d,'%6.3f'),'\pm',num2str(c1ed,'%6.3f'),'   offset=',num2str(c0d,'%6.2f'),'\pm',num2str(c0ed,'%6.2f')];
 nits=['Night fit  grad=',num2str(c1n,'%6.3f'),'\pm',num2str(c1en,'%6.3f'),'   offset=',num2str(c0n,'%6.2f'),'\pm',num2str(c0en,'%6.2f')];
 legend(h1,'Day (0900 - 1500 LT)','Night (2100 - 0300 LT)',days,nits,2)
 
 figure(2)
 subplot(2,2,1)
 h12=plot(tec(poid),latmd(poid),'.b',tec(poin),latmd(poin),'.r');
 grid 
 set(gca,'fontweight','bold')
 set(gca,'ytick',-90:30:90)
 xlabel('SAC-C TEC (TECU)','fontweight','bold')
 ylabel('Mag. Lat.','fontweight','bold')
 title('SAC-C','fontweight','bold')
 axis([0 ax3x -90 90])
 legend(h12,'Day','Night',4)
 
 subplot(2,2,2)
 h12=plot(tecgp(poid),latmd(poid),'.b',tecgp(poin),latmd(poin),'.r');
 grid 
 set(gca,'fontweight','bold')
 set(gca,'ytick',-90:30:90)
 xlabel('GPID TEC (TECU)','fontweight','bold')
 ylabel('Mag. Lat.','fontweight','bold')
 title(['GPID (',iris,')'],'fontweight','bold')
 axis([0 ax3x -90 90])
 legend(h12,'Day','Night',4)
 
 subplot(2,2,3)
 plot(tecgp,timmat,'.')
 grid
 set(gca,'ytick',0:4:24)
 set(gca,'fontweight','bold')
 xlabel('GPID TEC (TECU)','fontweight','bold')
 ylabel('Time (UT)','fontweight','bold')
 title(['UT (',iris,')'],'fontweight','bold')
 axis([0 ax1x 0 24])
 
 subplot(2,2,4)
 plot(tecgp,timl,'.')
 grid
 set(gca,'ytick',0:4:24)
 set(gca,'fontweight','bold')
 xlabel('GPID TEC (TECU)','fontweight','bold')
 ylabel('Time (LT)','fontweight','bold')
 title(['LT (SAC-C) (',iris,')'],'fontweight','bold')
 axis([0 ax1x 0 24])
 
elseif nopi==3
 
 elsd=rad2deg(els);

 dgv=10:10:80;
 nof=length(dgv);

 c1=nan*ones(1,nof);
 c0=c1;
 c1e=c1;
 c0e=c0;
 c1d=c1;
 c0d=c1;
 c1de=c1;
 c0de=c1;
 c1n=c1;
 c0n=c1;
 c1ne=c1;
 c0ne=c1;
 nopv=c1;
 nopvd=c1;
 nopvn=c1;
 lrg=c1;
 lrgd=c1;
 lrgn=c1;

 for i=1:nof

  poi=elsd>=dgv(i);
  poid=poi & (timl>=9 & timl<=15);
  poin=poi & (timl>=21 | timl<=3);
  
  nopv(i)=sum(poi);
  nopvd(i)=sum(poid);
  nopvn(i)=sum(poin);
  
  [c1(i),c0(i),c1e(i),c0e(i)]=linfit(tecgp(poi),tec(poi));
  [c1d(i),c0d(i),c1de(i),c0de(i)]=linfit(tecgp(poid),tec(poid));
  [c1n(i),c0n(i),c1ne(i),c0ne(i)]=linfit(tecgp(poin),tec(poin));
  
  lrg(i)=linrg(tecgp(poi),tec(poi));
  lrgd(i)=linrg(tecgp(poid),tec(poid));
  lrgn(i)=linrg(tecgp(poin),tec(poin));
  
 end

 figure(1)
 pwerr(dgv,c1,c1e,'g')
 hold on
 pwerr(dgv-0.5,c1d,c1de,'r')
 pwerr(dgv+0.5,c1n,c1ne,'b')
 h11=plot(dgv,c1,'gd',dgv-0.5,c1d,'rs',dgv+0.5,c1n,'bo');
 hold off
 set(h11(1),'markerfacecolor','g')
 set(h11(2),'markerfacecolor','r')
 set(h11(3),'markerfacecolor','b')
 set(gca,'fontweight','bold','fontsize',10)
 grid
 xlabel('SAC-C elevation mask angle (degrees)','fontweight','bold','fontsize',12)
 ylabel('Gradient','fontweight','bold','fontsize',12)
 tls=['Gradient: Various SAC-C mask (',iris,') Dno=',int2str(dno),' 1999'];
 title(tls,'fontweight','bold','fontsize',14)
 legend(h11,'Fit to full data','Day fit (0900-1500 LT)','Night fit (2100-0300 LT)',2)

 figure(2)
 pwerr(dgv,c0,c0e,'g')
 hold on
 pwerr(dgv,c0d,c0de,'r')
 pwerr(dgv,c0n,c0ne,'b')
 h12=plot(dgv,c0,'gd',dgv,c0d,'rs',dgv,c0n,'bo');
 hold off
 set(h12(1),'markerfacecolor','g')
 set(h12(2),'markerfacecolor','r')
 set(h12(3),'markerfacecolor','b')
 set(gca,'fontweight','bold','fontsize',10)
 grid
 xlabel('SAC-C elevation mask angle (degrees)','fontweight','bold','fontsize',12)
 ylabel('Offset (TECU)','fontweight','bold','fontsize',12)
 tls=['Offset: Various SAC-C mask (',iris,') Dno=',int2str(dno),' 1999'];
 title(tls,'fontweight','bold','fontsize',14)
 legend(h12,'Fit to full data','Day fit (0900-1500 LT)','Night fit (2100-0300 LT)',3)
 
 figure(3)
 subplot(2,1,1)
 h21=plot(dgv,nopv,'gd',dgv-0.5,nopvd,'rs',dgv+0.5,nopvn,'bo');
 set(h21(1),'markerfacecolor','g')
 set(h21(2),'markerfacecolor','r')
 set(h21(3),'markerfacecolor','b')
 set(gca,'fontweight','bold','fontsize',10)
 grid
 xlabel('SAC-C elevation mask angle (degrees)','fontweight','bold','fontsize',10)
 ylabel('Number of points','fontweight','bold','fontsize',10)
 tls=['Number of Points: Various SAC-C mask (',iris,') Dno=',int2str(dno),' 1999'];
 title(tls,'fontweight','bold','fontsize',12)
 legend(h21,'Fit to full data','Day fit (0900-1500 LT)','Night fit (2100-0300 LT)',1)
 
 subplot(2,1,2)
 h22=plot(dgv,lrg,'gd',dgv-0.5,lrgd,'rs',dgv+0.5,lrgn,'bo');
 set(h22(1),'markerfacecolor','g')
 set(h22(2),'markerfacecolor','r')
 set(h22(3),'markerfacecolor','b')
 set(gca,'fontweight','bold','fontsize',10)
 grid
 xlabel('SAC-C elevation mask angle (degrees)','fontweight','bold','fontsize',10)
 ylabel('Linear Regression','fontweight','bold','fontsize',10)
 tls=['Linear regression: Various SAC-C mask (',iris,') Dno=',int2str(dno),' 1999'];
 title(tls,'fontweight','bold','fontsize',12)
 legend(h22,'Fit to full data','Day fit (0900-1500 LT)','Night fit (2100-0300 LT)',1)
 
elseif nopi==4

 dgv=input('Minimum elevation angle to use from SAC-C (degrees)? ');

 elsd=rad2deg(els);

 latmb=0:10:90;
 nof=length(latmb)-1;

 c1=nan*ones(1,nof);
 c0=c1;
 c1e=c1;
 c0e=c1;
 c1d=c1;
 c0d=c1;
 c1de=c1;
 c0de=c1;
 c1n=c1;
 c0n=c1;
 c1ne=c1;
 c0ne=c1;
 nopv=c1;
 nopvd=c1;
 nopvn=c1;
 lrg=c1;
 lrgd=c1;
 lrgn=c1;
 
 for i=1:nof
  poim=abs(latmd)>=latmb(i) & abs(latmd)<latmb(i+1);
  poi=elsd>=dgv & poim;
  poid=poi & (timl>=9 & timl<=15);
  poin=poi & (timl>=21 | timl<=3);
  nopv(i)=sum(poi);
  nopvd(i)=sum(poid);
  nopvn(i)=sum(poin);
  if nopv(i)>2
   [c1(i),c0(i),c1e(i),c0e(i)]=linfit(tecgp(poi),tec(poi));
   lrg(i)=linrg(tecgp(poi),tec(poi));
  end
  if nopvd(i)>2
   [c1d(i),c0d(i),c1de(i),c0de(i)]=linfit(tecgp(poid),tec(poid));
   lrgd(i)=linrg(tecgp(poid),tec(poid));
  end
  if nopvn(i)>2
   [c1n(i),c0n(i),c1ne(i),c0ne(i)]=linfit(tecgp(poin),tec(poin));
   lrgn(i)=linrg(tecgp(poin),tec(poin));
  end
 end

 latmb2=5:10:85;

 figure(1)
 pwerr(latmb2,c1,c1e,'g')
 hold on
 pwerr(latmb2-0.5,c1d,c1de,'r')
 pwerr(latmb2+0.5,c1n,c1ne,'b')
 h11=plot(latmb2,c1,'gd',latmb2-0.5,c1d,'rs',latmb2+0.5,c1n,'bo');
 hold off
 axis([0 90 0 2])
 set(h11(1),'markerfacecolor','g')
 set(h11(2),'markerfacecolor','r')
 set(h11(3),'markerfacecolor','b')
 set(gca,'fontweight','bold','fontsize',10)
 grid
 xlabel('SAC-C magnetic latitude (degrees)','fontweight','bold','fontsize',12)
 ylabel('Gradient','fontweight','bold','fontsize',12)
 tls=['Gradient: ',int2str(dgv),'^o SAC-C mask (',iris,') Dno=',int2str(dno),' 1999'];
 title(tls,'fontweight','bold','fontsize',14)
 legend(h11,'Fit to full data','Day fit (0900-1500 LT)','Night fit (2100-0300 LT)',2)
 for j=1:9
  ttv=cell(3,1);
  ttv(1)={['F ',num2str(round(100*lrg(j))/100),'  ',int2str(nopv(j))]};
  ttv(2)={['D ',num2str(round(100*lrgd(j))/100),'  ',int2str(nopvd(j))]};
  ttv(3)={['N ',num2str(round(100*lrgn(j))/100),'  ',int2str(nopvn(j))]};
  text(latmb2(j)-3.5,0.15,ttv,'fontweight','bold','color','r','fontsize',8)
 end
 
 figure(2)
 pwerr(latmb2,c0,c0e,'g')
 hold on
 pwerr(latmb2,c0d,c0de,'r')
 pwerr(latmb2,c0n,c0ne,'b')
 h12=plot(latmb2,c0,'gd',latmb2,c0d,'rs',latmb2,c0n,'bo');
 hold off
 set(h12(1),'markerfacecolor','g')
 set(h12(2),'markerfacecolor','r')
 set(h12(3),'markerfacecolor','b')
 set(gca,'fontweight','bold','fontsize',10)
 grid
 xlabel('SAC-C magnetic latitude (degrees)','fontweight','bold','fontsize',12)
 ylabel('Offset (TECU)','fontweight','bold','fontsize',12)
 tls=['Offset: ',int2str(dgv),'^o SAC-C mask (',iris,') Dno=',int2str(dno),' 1999'];
 title(tls,'fontweight','bold','fontsize',14)
 legend(h12,'Fit to full data','Day fit (0900-1500 LT)','Night fit (2100-0300 LT)',2)
 
else
 
 poid=timl>=9 & timl<=15;
 poin=timl>=21 | timl<=3;
 
 figure(1)
 dtec=tec-tecgp;
 plot(timl,dtec,'b.')
 grid
 tecl=2*floor(min(dtec)/2);
 tecu=2*ceil(max(dtec)/2);
 set(gca,'xtick',0:2:24)
 set(gca,'ytick',tecl:2:tecu)
 axis([0 24 tecl tecu])
 set(gca,'fontweight','bold','fontsize',12)
 xlabel('SAC-C Local Time','fontweight','bold','fontsize',12)
 ylabel('GPID SAC-C TEC (TECU)','fontweight','bold','fontsize',12)
 tls=['TEC Difference (',iris,') Dno=',int2str(dno),' 1999'];
 title(tls,'fontweight','bold','fontsize',14)
 
 figure(2)
 h2=plot(latmd(poid),dtec(poid),'b.',latmd(poin),dtec(poin),'r.');
 grid
 tecl=2*floor(min(dtec)/2);
 tecu=2*ceil(max(dtec)/2);
 set(gca,'xtick',-90:15:90)
 set(gca,'ytick',tecl:2:tecu)
 axis([-90 90 tecl tecu])
 set(gca,'fontweight','bold','fontsize',12)
 xlabel('SAC-C Magnetic Latitude','fontweight','bold','fontsize',12)
 ylabel('GPID SAC-C TEC (TECU)','fontweight','bold','fontsize',12)
 tls=['TEC Difference (',iris,') Dno=',int2str(dno),' 1999'];
 title(tls,'fontweight','bold','fontsize',14)
 legend(h2,'Day (0900 - 1500 LT)','Night (2100 - 0300 LT)',3)

end