%% pwtpfit8.m is an updated version of the ealier fliptp6.m.
%% Its purpose is to fit Titheridge's temperature function to a variety of FLIP
%% temperature profiles, so that the variation in T0 and G0 as function of F10.7
%% can be estimated. For a given L-shell, these variations are assumed to be
%% linear with F10.7. In each case the linear coefficients and their errors are 
%% estimated. Using the errors as weightings, the coefficients are then fitted
%% with polynominals to determine their variation with latitude.
%%

clear all
close all

lonlsh=deg2rad(70);
lshv=[1.25 1.5 2 3 4 6 8];
nop=length(lshv);

%Result vectors for co-efficients of initial fits (value and error)
ptd1=zeros(nop,2);            %day base temperature
ptd0=ptd1;
ptn1=ptd1;                    %night base temperature
ptn0=ptd1;
pgd1=ptd1;                    %day base gradient
pgd0=ptd1;
pgn1=ptd1;                    %night base gradient
pgn0=ptd1;
t0dm=[];
t0nm=[];
g0dm=[];
g0nm=[];
t0dsm=[];
t0nsm=[];
g0dsm=[];
g0nsm=[];

%Loop through each L-shell, determining co-efficents of required fits
if exist('pwtpfit8.mat')==2
 load pwtpfit8
else
 tic
 for i=1:nop
  [i nop]

  %Fit Titheridge equation to FLIP results
  [t0d,g0d,t0n,g0n,t0ds,g0ds,t0ns,g0ns,edd,edn,kpm,f107m,dnov,yearv]=fliptpf4(lonlsh,lshv(i));

  %Conduct linear fit to T0 and G0 values
  [ptd1i,ptd0i,ptd1ei,ptd0ei]=linfit2(f107m,t0d,t0ds);
  [ptn1i,ptn0i,ptn1ei,ptn0ei]=linfit2(f107m,t0n,t0ns);
  [pgd1i,pgd0i,pgd1ei,pgd0ei]=linfit2(f107m,g0d,g0ds);
  [pgn1i,pgn0i,pgn1ei,pgn0ei]=linfit2(f107m,g0n,g0ns);
 
  %Record results of fit
  ptd1(i,:)=[ptd1i ptd1ei];    
  ptd0(i,:)=[ptd0i ptd0ei];
  ptn1(i,:)=[ptn1i ptn1ei];
  ptn0(i,:)=[ptn0i ptn0ei];
  pgd1(i,:)=[pgd1i pgd1ei];
  pgd0(i,:)=[pgd0i pgd0ei];
  pgn1(i,:)=[pgn1i pgn1ei];
  pgn0(i,:)=[pgn0i pgn0ei];
 
  %Record the values fits conducted to
  t0dm=[t0dm;t0d];
  t0nm=[t0nm;t0n];
  g0dm=[g0dm;g0d];
  g0nm=[g0nm;g0n];
  t0dsm=[t0dsm;t0ds];
  t0nsm=[t0nsm;t0ns];
  g0dsm=[g0dsm;g0ds];
  g0nsm=[g0nsm;g0ns];

  toc
 
 end 
 
 save pwtpfit8
 
end

re=6371;
sv=ones(1,nop);
latm=lshpar(lshv,-sv,re+350*sv);
sl=sin(latm').^2;
slf=linspace(0,1,100);

%Conduct fits to fitted linear coefficients using a weighted polynomial fit
pptd1=pwpolfit(sl,ptd1(:,1),ptd1(:,2),3);         
pptd0=pwpolfit(sl,ptd0(:,1),ptd0(:,2),3);         
pptn1=pwpolfit(sl,ptn1(:,1),ptn1(:,2),3);         
pptn0=pwpolfit(sl,ptn0(:,1),ptn0(:,2),3);         
ppgd1=pwpolfit(sl,pgd1(:,1),pgd1(:,2),3);         
ppgd0=pwpolfit(sl,pgd0(:,1),pgd0(:,2),3);         
ppgn1=pwpolfit(sl,pgn1(:,1),pgn1(:,2),3);         
ppgn0=pwpolfit(sl,pgn0(:,1),pgn0(:,2),3);         

pptd1f=polyval(pptd1,slf);
pptd0f=polyval(pptd0,slf);
pptn1f=polyval(pptn1,slf);
pptn0f=polyval(pptn0,slf);
ppgd1f=polyval(ppgd1,slf);
ppgd0f=polyval(ppgd0,slf);
ppgn1f=polyval(ppgn1,slf);
ppgn0f=polyval(ppgn0,slf);

save pwf107 pptd1 pptd0 pptn1 pptn0 ppgd1 ppgd0 ppgn1 ppgn0

disp(' ')
fprintf('pptd1=[%6.3e %6.3e %6.3e %6.3e]\n',pptd1)
fprintf('pptd0=[%6.3e %6.3e %6.3e %6.3e]\n',pptd0)
fprintf('pptn1=[%6.3e %6.3e %6.3e %6.3e]\n',pptn1)
fprintf('pptn0=[%6.3e %6.3e %6.3e %6.3e]\n',pptn0)
disp(' ')
fprintf('ppgd1=[%6.3e %6.3e %6.3e %6.3e]\n',ppgd1)
fprintf('ppgd0=[%6.3e %6.3e %6.3e %6.3e]\n',ppgd0)
fprintf('ppgn1=[%6.3e %6.3e %6.3e %6.3e]\n',ppgn1)
fprintf('ppgn0=[%6.3e %6.3e %6.3e %6.3e]\n',ppgn0)
disp(' ')

f107f=linspace(60,250);
pli=0;

for k=1:4
 figure(k)
 pli=pli+1;
 subplot(2,2,1)
 h1=plot(f107m,t0dm(pli,:),'ob',f107f,ptd1(pli,1)*f107f+ptd0(pli,1),'b',f107m,t0nm(pli,:),'^r',f107f,ptn1(pli,1)*f107f+ptn0(pli,1),'r');
 hold on
 pwerr(f107m,t0dm(pli,:),t0dsm(pli,:),'b')
 pwerr(f107m,t0nm(pli,:),t0nsm(pli,:),'r')
 hold off
 axis([50 250 500 3500])
 set(h1(1),'markerfacecolor','b')
 set(h1([2 3 4]),'linewidth',1.5)
 set(gca,'fontweight','bold','fontsize',14)
 grid
 xlabel('F_{10.7}','fontweight','bold','fontsize',14)
 ylabel('{\itT_o} (K)','fontweight','bold','fontsize',14)
 title(['L=',num2str(lshv(pli))],'fontweight','bold','fontsize',14)
 legend(h1([1 3]),'Day','Night',1)

 subplot(2,2,3)
 h1=plot(f107m,g0dm(pli,:),'ob',f107f,pgd1(pli,1)*f107f+pgd0(pli,1),'b',f107m,g0nm(pli,:),'^r',f107f,pgn1(pli,1)*f107f+pgn0(pli,1),'r');
 hold on
 pwerr(f107m,g0dm(pli,:),g0dsm(pli,:),'b')
 pwerr(f107m,g0nm(pli,:),g0nsm(pli,:),'r')
 hold off
 axis([50 250 -1 8])
 set(h1(1),'markerfacecolor','b')
 set(h1([2 3 4]),'linewidth',1.5)
 set(gca,'fontweight','bold','fontsize',14)
 grid
 xlabel('F_{10.7}','fontweight','bold','fontsize',14)
 ylabel('{\itG_o} (K/km)','fontweight','bold','fontsize',14)
 %title(['L=',num2str(lshv(pli))],'fontweight','bold','fontsize',12)

 if k<4
  pli=pli+1;
  subplot(2,2,2)
  h1=plot(f107m,t0dm(pli,:),'ob',f107f,ptd1(pli,1)*f107f+ptd0(pli,1),'b',f107m,t0nm(pli,:),'^r',f107f,ptn1(pli,1)*f107f+ptn0(pli,1),'r');
  hold on
  pwerr(f107m,t0dm(pli,:),t0dsm(pli,:),'b')
  pwerr(f107m,t0nm(pli,:),t0nsm(pli,:),'r')
  hold off
  axis([50 250 500 3500])
  set(h1(1),'markerfacecolor','b')
  set(h1([2 3 4]),'linewidth',1.5)
  set(gca,'fontweight','bold','fontsize',14)
  grid
  xlabel('F_{10.7}','fontweight','bold','fontsize',14)
  ylabel('{\itT_o} (K)','fontweight','bold','fontsize',14)
  title(['L=',num2str(lshv(pli))],'fontweight','bold','fontsize',14)

  subplot(2,2,4)
  h1=plot(f107m,g0dm(pli,:),'ob',f107f,pgd1(pli,1)*f107f+pgd0(pli,1),'b',f107m,g0nm(pli,:),'^r',f107f,pgn1(pli,1)*f107f+pgn0(pli,1),'r');
  hold on
  pwerr(f107m,g0dm(pli,:),g0dsm(pli,:),'b')
  pwerr(f107m,g0nm(pli,:),g0nsm(pli,:),'r')
  hold off
  axis([50 250 -1 8])
  set(h1(1),'markerfacecolor','b')
  set(h1([2 3 4]),'linewidth',1.5)
  set(gca,'fontweight','bold','fontsize',14)
  grid
  xlabel('F_{10.7}','fontweight','bold','fontsize',14)
  ylabel('{\itG_o} (K/km)','fontweight','bold','fontsize',14)
  %title(['L=',num2str(lshv(pli))],'fontweight','bold','fontsize',12)
 end
end

figure(5)
subplot(2,2,1)
h11=plot(sl,ptd1(:,1),'ob',slf,pptd1f,'b',sl+0.005,ptn1(:,1),'r^',slf,pptn1f,'r');
hold on
pwerr(sl,ptd1(:,1),ptd1(:,2),'b')
pwerr(sl+0.005,ptn1(:,1),ptn1(:,2),'r')
hold off
grid
set(h11(1),'markerfacecolor','b')
set(h11([2 3 4]),'linewidth',1.5)
set(gca,'xtick',0:0.2:1)
xlabel('{\its_L}','fontweight','bold','fontsize',14)
ylabel('{\itT}_{01} (K/F_{10.7})','fontweight','bold','fontsize',14)
set(gca,'fontweight','bold','fontsize',14)

subplot(2,2,2)
h11=plot(sl,ptd0(:,1),'ob',slf,pptd0f,'b',sl+0.005,ptn0(:,1),'r^',slf,pptn0f,'r');
hold on
pwerr(sl,ptd0(:,1),ptd0(:,2),'b')
pwerr(sl+0.005,ptn0(:,1),ptn0(:,2),'r')
hold off
grid
set(h11(1),'markerfacecolor','b')
set(h11([2 3 4]),'linewidth',1.5)
set(gca,'xtick',0:0.2:1)
xlabel('{\its_L}','fontweight','bold','fontsize',14)
ylabel('{\itT}_{00} (K)','fontweight','bold','fontsize',14)
set(gca,'fontweight','bold','fontsize',14)

subplot(2,2,3)
h11=plot(sl,pgd1(:,1),'ob',slf,ppgd1f,'b',sl+0.005,pgn1(:,1),'r^',slf,ppgn1f,'r');
hold on
pwerr(sl,pgd1(:,1),pgd1(:,2),'b')
pwerr(sl+0.005,pgn1(:,1),pgn1(:,2),'r')
hold off
grid
set(h11(1),'markerfacecolor','b')
set(h11([2 3 4]),'linewidth',1.5)
set(gca,'xtick',0:0.2:1)
xlabel('{\its_L}','fontweight','bold','fontsize',14)
ylabel('{\itG}_{01} (K/km/F_{10.7})','fontweight','bold','fontsize',14)
set(gca,'fontweight','bold','fontsize',14)

subplot(2,2,4)
h11=plot(sl,pgd0(:,1),'ob',slf,ppgd0f,'b',sl+0.005,pgn0(:,1),'r^',slf,ppgn0f,'r');
hold on
pwerr(sl,pgd0(:,1),pgd0(:,2),'b')
pwerr(sl+0.005,pgn0(:,1),pgn0(:,2),'r')
hold off
grid
set(h11(1),'markerfacecolor','b')
set(h11([2 3 4]),'linewidth',1.5)
set(gca,'xtick',0:0.2:1)
xlabel('{\its_L}','fontweight','bold','fontsize',14)
ylabel('{\itG}_{00} (K/km)','fontweight','bold','fontsize',14)
set(gca,'fontweight','bold','fontsize',14)
legend(h11([1 3]),'Day','Night',4)
