function out = son_comb75z
%%
%% This program is written for combining 75deg elev and zenith Sonde ISR data in ascii
%% files made from rd_son.m
%% BA Emery, April 2011
clear all; close all
warning off MATLAB:divideByZero

	glat = 66.99;
	glon = -50.95;

%  kindat=5501 long-pulse ~4.5min 3-position (Z, NE+NW 75elev) 98-1063/95-1053km
%    Linearly interpolate heights from 2 75deg elev 20-min medians from 22 hts 95-1053 km
%     to 21 hts of Z profile from 98-1063 km.

  i20w2 = sprintf ('son_lp050514_15_20m_profiles');
  INI20W2 = fopen (i20w2,'r');
  k20w2 = sprintf ('son_lp75_050514_15_20m_profiles');
  INK20W2 = fopen (k20w2,'r');
  j20w2 = sprintf ('son_lpc_050514_15_20m_profiles');
  j20w1 = sprintf ('son_lpc_050514_15_20m_nmax');
  OUTK20W2 = fopen (j20w2,'w');
 fprintf (OUTK20W2,' year dayn   uth alt(km) #Ti medTi(K)    stdTi #Te medTe(K)    stdTe #Ne   medNe(m-3)   stdNe \n');
  OUTK20W1 = fopen (j20w1,'w');
 fprintf (OUTK20W1,' year dayn moda   uth   glat   glon   hmax foF2MHz  Nmax(m-3)    TEC(m-2) \n');

% Set altitude and Ne and fill in gaps in Ne (jne) above 100 km (j100)
% Values for kindat=5501
	j100 = 1;
	j100m1 = min(1,j100-1);
	j100p = max(j100,2);
% Look at peak >100km instead of >150km
	j150 = j100;
% Delete 98km for long-pulse (or 95 km for 75 deg elev)
	j150 = 2;
% Consider hmax only below 600 km; LP htz(13) = 560km, htjm(14)=563km for 75 deg elev
j600 = 13;
j600m1 = j600 -1;

ialt = 21;
kalt = 22;
n20m = 48;
missing = -32767.
val1 = 1.;

% Skip over first line
 linobsi = fgetl(INI20W2);
 linobsk = fgetl(INK20W2);


nti = 6;
nte = 9;
nne = 12;
for m=1:n20m
% Read 2D profiles
 for n9=1:ialt
   linobsi = fgetl(INI20W2);
   obsidat = sscanf(linobsi,'%f %f %f %f %f %f %f %f %f %f %f %f %f');
if (n9==1)
 ibyrt = obsidat(1);
 dayi(m) = obsidat(2);
 uthi(m) = obsidat(3);
end
   alti(n9) = obsidat(4);
   htjm(n9) = obsidat(4);
   numtii(n9,m) = obsidat(nti-1);
   tii(n9,m) = obsidat(nti);
   stdtii(n9,m) = obsidat(nti+1);
   numtei(n9,m) = obsidat(nte-1);
   tei(n9,m) = obsidat(nte);
   stdtei(n9,m) = obsidat(nte+1);
   numnei(n9,m) = obsidat(nne-1);
   nei(n9,m) = obsidat(nne);
   stdnei(n9,m) = obsidat(nne+1);
 end	% end n9=1:ialt
 for n9=1:kalt
   linobsk = fgetl(INK20W2);
   obsidat = sscanf(linobsk,'%f %f %f %f %f %f %f %f %f %f %f %f %f');
if (n9==1)
 dayk(m) = obsidat(2);
 uthk(m) = obsidat(3);
end
   altk(n9) = obsidat(4);
   numtik(n9,m) = obsidat(nti-1);
   tik(n9,m) = obsidat(nti);
   stdtik(n9,m) = obsidat(nti+1);
   numtek(n9,m) = obsidat(nte-1);
   tek(n9,m) = obsidat(nte);
   stdtek(n9,m) = obsidat(nte+1);
   numnek(n9,m) = obsidat(nne-1);
   nek(n9,m) = obsidat(nne);
   stdnek(n9,m) = obsidat(nne+1);
 end	% end n9=1:kalt

if (m==1)
% Linearly interpolate in height, except highest altitude #22k 1053 give to #21z 1068 km
% calculate linear interpolations between 2 altitude regions #15z=#16k at 686km
% Use fac1 for j, and fac2 for j+1
% #1-14z (99-623km)
 for j=1:14
   fac1(j) = (altk(j+1)-alti(j))/(altk(j+1)-altk(j));
   fac2(j) = (alti(j)-altk(j))/(altk(j+1)-altk(j));
 end
fac1(15) = 1.;
fac2(15) = 1.;
fac1(21) = 1.;
fac2(21) = 1.;
% #16-20z (749-1001km)
 for j=16,kalt-2
   fac1(j) = (altk(j+2)-alti(j))/(altk(j+2)-altk(j+1));
   fac2(j) = (alti(j)-altk(j+1))/(altk(j+2)-altk(j+1));
 end
end	% if m==1

avut = 0.5*(uthi(m) + uthk(m));
mdayn = dayi(m);
moda = 514;
if (mdayn== 135)
 moda = 515;
end
if (dayi(m)<dayk(m))
  avut = 0.5*(uthi(m)+uthk(m)+24.);
end

for j=1:14
% Calculate average values weighted by numbers if have more than 1 in bin
 [numtij,avtij,stdtij] = htinterp(numtik(j,m),numtik(j+1,m),tik(j,m),tik(j+1,m),stdtik(j,m),stdtik(j+1,m),fac1(j),fac2(j),missing);
 [numti,avti(j),averti(j)] = htinterp(numtii(j,m),numtij,tii(j,m),avtij,stdtii(j,m),stdtij,val1,val1,missing);
 [numtij,avtij,stdtij] = htinterp(numtek(j,m),numtek(j+1,m),tek(j,m),tek(j+1,m),stdtek(j,m),stdtek(j+1,m),fac1(j),fac2(j),missing);
 [numte,avte(j),averte(j)] = htinterp(numtei(j,m),numtij,tei(j,m),avtij,stdtei(j,m),stdtij,val1,val1,missing);
 [numtij,avtij,stdtij] = htinterp(numnek(j,m),numnek(j+1,m),nek(j,m),nek(j+1,m),stdnek(j,m),stdnek(j+1,m),fac1(j),fac2(j),missing);
 [numne,avne(j),averne(j)] = htinterp(numnei(j,m),numtij,nei(j,m),avtij,stdnei(j,m),stdtij,val1,val1,missing);
fprintf (OUTK20W2,' %4d %4d %5.2f %7.2f %3d %8.1f %8.1f %3d %8.1f %8.1f %3d %11.4e %11.4e \n',ibyrt,mdayn,avut,alti(j),numti,avti(j),averti(j),numte,avte(j),averte(j),numne,avne(j),averne(j) ); 
end
j=15;
 [numti,avti(j),averti(j)] = htinterp(numtii(j,m),numtik(j+1,m),tii(j,m),tik(j+1,m),stdtii(j,m),stdtik(j+1,m),val1,val1,missing);
 [numte,avte(j),averte(j)] = htinterp(numtei(j,m),numtek(j+1,m),tei(j,m),tek(j+1,m),stdtei(j,m),stdtek(j+1,m),val1,val1,missing);
 [numne,avne(j),averne(j)] = htinterp(numnei(j,m),numnek(j+1,m),nei(j,m),nek(j+1,m),stdnei(j,m),stdnek(j+1,m),val1,val1,missing);
fprintf (OUTK20W2,' %4d %4d %5.2f %7.2f %3d %8.1f %8.1f %3d %8.1f %8.1f %3d %11.4e %11.4e \n',ibyrt,mdayn,avut,alti(j),numti,avti(j),averti(j),numte,avte(j),averte(j),numne,avne(j),averne(j) ); 
for j=16:20
 [numtij,avtij,stdtij] = htinterp(numtik(j,m),numtik(j+2,m),tik(j,m),tik(j+2,m),stdtik(j,m),stdtik(j+2,m),fac1(j),fac2(j),missing);
 [numti,avti(j),averti(j)] = htinterp(numtii(j,m),numtij,tii(j,m),avtij,stdtii(j,m),stdtij,val1,val1,missing);
 [numtij,avtij,stdtij] = htinterp(numtek(j,m),numtek(j+2,m),tek(j,m),tek(j+2,m),stdtek(j,m),stdtek(j+2,m),fac1(j),fac2(j),missing);
 [numte,avte(j),averte(j)] = htinterp(numtei(j,m),numtij,tei(j,m),avtij,stdtei(j,m),stdtij,val1,val1,missing);
 [numtij,avtij,stdtij] = htinterp(numnek(j,m),numnek(j+2,m),nek(j,m),nek(j+2,m),stdnek(j,m),stdnek(j+2,m),fac1(j),fac2(j),missing);
 [numne,avne(j),averne(j)] = htinterp(numnei(j,m),numtij,nei(j,m),avtij,stdnei(j,m),stdtij,val1,val1,missing);
fprintf (OUTK20W2,' %4d %4d %5.2f %7.2f %3d %8.1f %8.1f %3d %8.1f %8.1f %3d %11.4e %11.4e \n',ibyrt,mdayn,avut,alti(j),numti,avti(j),averti(j),numte,avte(j),averte(j),numne,avne(j),averne(j) ); 
end
j=21;
 [numti,avti(j),averti(j)] = htinterp(numtii(j,m),numtik(j+1,m),tii(j,m),tik(j+1,m),stdtii(j,m),stdtik(j+1,m),val1,val1,missing);
 [numte,avte(j),averte(j)] = htinterp(numtei(j,m),numtek(j+1,m),tei(j,m),tek(j+1,m),stdtei(j,m),stdtek(j+1,m),val1,val1,missing);
 [numne,avne(j),averne(j)] = htinterp(numnei(j,m),numnek(j+1,m),nei(j,m),nek(j+1,m),stdnei(j,m),stdnek(j+1,m),val1,val1,missing);
fprintf (OUTK20W2,' %4d %4d %5.2f %7.2f %3d %8.1f %8.1f %3d %8.1f %8.1f %3d %11.4e %11.4e \n',ibyrt,mdayn,avut,alti(j),numti,avti(j),averti(j),numte,avte(j),averte(j),numne,avne(j),averne(j) ); 

%  Calculate hmF2, NmF2, and TEC (1 TECU=10^16elec/m2)
	xnemax = 0.;
	for j=j150:j600m1
%  Find max Ne 
	  if (avne(j) > xnemax)
	    xnemax = avne(j);
	    jx = j;
	  end
        end
	h(2) = htjm(jx);
	p(2) = avne(jx);
	hmf2 = htjm(jx);
	xnmax = avne(jx);
%  Check to see if max Ne next to missing data or not
	for j=jx-1:-1:j150
	  if (avne(j) > 0.), break, end
        end
       if (j > j150)
	h(1) = htjm(j);
	p(1) = avne(j);
	for j=jx+1:j600
	  if (avne(j) > 0.), break, end
        end
       if (j < j600)
	h(3) = htjm(j);
	p(3) = avne(j);
%  Use polynomial to find hmF2 and NmF2
        [hmf2,xnmax] = hnmf2 (h,p);
       end
       end
      fof2 = sqrt(xnmax/1.24e+10);

%  Find tec from non-missing data
	tec = 0.;
	htjm1 = htjm(max(1,j100m1));
        for j=j100p:ialt
	  if (avne(j) > 0.) 
	    tec = tec + avne(j)*(htjm(j)-htjm1)*1.e+3;
	    htjm1 = htjm(j);
	  end
        end

fprintf (OUTK20W1,' %4d %4d %4d %5.2f %6.2f %6.2f %7.2f %5.2f %11.4e %11.4e \n',ibyrt, mdayn,moda,avut,glat,glon,hmf2,fof2,xnmax,tec);

end	% end of for m=1:n20m

function [num3,av3,std3] = htinterp(num1,num2,av1,av2,std1,std2,fac1,fac2,missing)
if (num1 > 1 & num2 > 1)
 num3 = fac1*num1 + fac2*num2;
 av3 = (fac1*num1*av1 + fac2*num2*av2)/num3;
 std3 = (fac1*num1*std1 + fac2*num2*std2)/num3;
elseif (num1 > 1 & num2 <=1 & fac1 >= fac2)
 num3 = num1;
 av3 = av1;
 std3 = std1;
elseif (num1 <=1 & num2 > 1 & fac2 >= fac1)
 num3 = num2;
 av3 = av2;
 std3 = std2;
else
 num3 = 0;
 av3 = missing;
 std3 = missing;
end

return

%ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
function [hmf2,nmf2] = hnmf2(h,p)

% Find the Nmax (F2 or E) and hmax based on the electron profile.
% btf: This is taken from an idl procedure provided by Liying Qian 
% bae 4/11: Ben Foster's Fortran version was revised to matlab 

        h12 = h(1)*h(1);
        h22 = h(2)*h(2);
        h32 = h(3)*h(3);
        deltx=h12*h(2)+h22*h(3)+h32*h(1)-h32*h(2)-h12*h(3)-h22*h(1);
        atx=p(1)*h(2)+p(2)*h(3)+p(3)*h(1)-h(2)*p(3)-h(3)*p(1)-h(1)*p(2);
        ax=atx/deltx;
        btx=h12*p(2)+h22*p(3)+h32*p(1)-h32*p(2)-h12*p(3)-h22*p(1);
        bx=btx/deltx;
        ctx=h12*h(2)*p(3)+h22*h(3)*p(1)+h32*h(1)*p(2)-h32*h(2)*p(1)-h12*h(3)*p(2)-h22*h(1)*p(3);
        cx=ctx/deltx;
        hmf2=-(bx/(2.*ax));
        nmf2=-((bx*bx-4.*ax*cx)/(4.*ax));
 
return
