
function out = rd_eson
%%
%% This program is written for reading Sonde ISR ExB CEDAR DB .asc data
%% BA Emery, May 2011
clear all; close all
warning off MATLAB:divideByZero
 
%  kindat=15011 EPEC (not 3km) from long-pulse ~4.5min 3-position (Z, NE+NW 75elev) every 1.3min over ~4min but saved every 20 min
%  Use EperpE,N 1640,1650 with |B|=0.4884 gauss at 298km to calc ViperpN,E, then Ui,Vi,Wi
%  Save also averaged 176-563km ViperpE,N if error bar<100m/s for cf
%  Program to read the ascii format data submitted to the CEDAR DB
%   .asc not available from DB access system directly, but can be converted from .cbf

	missing = -32767;

% Open 1 file for 1D |B|298km, EperpE,N, ViperE,N298km, Ui,Vi,Wi (+/-40km?)

% EPEC (not 3 km) 
% Set previous minute to same as first time read
	k20 = 1;
	mp(1) = 3;
	mp(2) = 3;
%  datfil = sprintf ('/home/emery/archive/son050514g.010_epec');
	mp(1) = 2;
	mp(2) = 2;
  datfil = sprintf ('/home/emery/archive/son061214g.014_epec_3p');
  DAT = fopen (datfil,'r');

	nm(1) = 0;
	nm(2) = 0;
	glat = 66.99;
	glon = -50.95;
	bmag298 = 0.4884
% 1-D EperpE,N and errs at 9,11 and 10,12
	iepe = 9;
	ierepe = 10;
	iepn = 11;
	ierepn = 12;
% 2-D Bmag; ViperpE,N,Par and errs at 3; 4,6,8 and 5,7,9
	jb = 3;
	jvpe = 4;
	jervpe = 5;
	jvpn = 6;
	jervpn = 7;
	jvpp = 8;
	jervpp = 9;
% Look at aves >166km
	j150 = 4;

%  iw1 = sprintf ('son_epec_050514_15_ExB');
  iw1 = sprintf ('son_epec_061214_15_ExB');
  OUTIW1 = fopen (iw1,'w');
 fprintf (OUTIW1,' year dayn moda   uth   glat   glon  ViperpE,N,par m/s (176-563km)  Bgauss_298km EperpE,err,ViExBn  EperpN,err,ViExBe \n');
%  i20w1 = sprintf ('son_epec_050514_15_20m_ExB');
  i20w1 = sprintf ('son_epec_061214_15_20m_ExB');
  OUTI20W1 = fopen (i20w1,'w');
 fprintf (OUTI20W1,' year dayn moda   uth   glat   glon  #,val,err_ViperpE,N,par m/s (176-563km)  Bgauss_298km #,EperpE,err,ViExBn  #,EperpN,err,ViExBe \n');


        nrecrd = 0;

while(feof(DAT)==0)
     
% Read Prologue
        nrecrd = nrecrd + 1;
	line = fgetl(DAT);
	lprol = sscanf(line,'%d');

%  Read (and write) text portions (headers 3101 or catalogs 2101)
 if (lprol(2)==3101 | lprol(2)==2101)
  if (nrecrd==1) 
   lprol
  end
      for n=2:lprol(1)
	line = fgetl(DAT);
	char80 = sscanf(line,'%s');
      end
 else
        kinst = lprol(3);
        kindat = lprol(4);
        ibyrt = lprol(5);
        ibdtt = lprol(6);
        ibhmt = lprol(7);
        ibcst = lprol(8);
        ieyrt = lprol(9);
        iedtt = lprol(10);
        iehmt = lprol(11);
        iecst = lprol(12);
        jpar = lprol(14);
        mpar = lprol(15);
        nrows = lprol(16);

% Convert year mmdd hhmm csec to daynumber ndayn
        itime(1) = ibyrt;
        itime(2) = ibdtt;
        itime(3) = ibhmt;
        itime(4) = ibcst;
        [ndayn] = cvt2dn (itime);

%  Read JPAR 1-D parameters, and NROWS of MPAR 2-D parameters from instrument KINST
%   between begin date IBYRT,IBDTT,IBHMT,IBCST and end date IEYRT,IEDTT,IEHMT,IECST
%   (year, monthday, UThourmin, centisecond)
%      write (6,"(1x,'ndayn kinst kindat begin,end_time 1-D 2-D rows=',
%    |   14i6)") ndayn, (lprol(i1),i1=3,12),(lprol(i2),i2=14,16)
        if (jpar > 60 | mpar > 60 | nrows > 900) 
'large ',jpar,mpar,nrows
% The stop will make it stop abruptly with an error to fix jpar,mpar>60 etc
stop
        end
% Are there any JPAR (single valued parameters)?
       if (jpar > 0) 
%  Read 1-D parameter codes
	line = fgetl(DAT);
	i20a = sscanf(line,'%d');
	for i=1:min(jpar,20)
	  icod1d(i) = i20a(i);
	end
	if (jpar > 20)
	 i20b = sscanf(line,'%d');
	 for i=21:min(jpar,40)
	  icod1d(i) = i20b(i-20);
	 end
	 if (jpar > 40)
	  i20c = sscanf(line,'%d');
	  for i=41:min(jpar,60)
	   icod1d(i) = i20c(i-40);
	  end
	 end
	end	  
%  Read 1-D parameter values
	line = fgetl(DAT);
	i20a = sscanf(line,'%d');
	for i=1:min(jpar,20)
	  ipar1d(i) = i20a(i);
	end
	if (jpar > 20)
	 i20b = sscanf(line,'%d');
	 for i=21:min(jpar,40)
	  ipar1d(i) = i20b(i-20);
	 end
	end	
	 if (jpar > 40)
	  i20c = sscanf(line,'%d');
	  for i=41:min(jpar,60)
	   ipar1d(i) = i20c(i-40);
	  end
	 end
       end

% Are there any MPAR (multiple valued parameters)?
   if (mpar > 0) 
%  Read 2-D parameter codes
	line = fgetl(DAT);
	j20a = sscanf(line,'%d');
	for i=1:min(mpar,20)
	  icod2d(i) = j20a(i);
	end
	if (mpar > 20)
	 line = fgetl(DAT);
	 j20b = sscanf(line,'%d');
	 for i=21:min(mpar,40)
	  icod2d(i) = j20b(i-20);
	 end
	 if (mpar > 40)
	  line = fgetl(DAT);
	  j20c = sscanf(line,'%d');
	  for i=41:min(mpar,60)
	   icod2d(i) = j20c(i-40);
	  end
	 end
	end	  
%  Read 2-D parameter values
clear ipar2d;
  for j=1:nrows
	line = fgetl(DAT);
	j20a = sscanf(line,'%d');
	for i=1:min(mpar,20)
	  ipar2d(i,j) = j20a(i);
	end
	if (mpar > 20)
	 line = fgetl(DAT);
	 j20b = sscanf(line,'%d');
	 for i=21:min(mpar,40)
	  ipar2d(i,j) = j20b(i-20);
	 end
	 if (mpar > 40)
	  line = fgetl(DAT);
	  j20c = sscanf(line,'%d');
	  for i=41:min(mpar,60)
	   ipar2d(i,j) = j20c(i-40);
	  end
	 end
	end	  
  end


	ihr = fix(ibhmt/100);
	imn = fix(ibhmt - ihr*100);
	uth = ihr + imn/60.;
	glat = ipar1d(1)*0.01;
	glon = ipar1d(2)*0.01;

% Calculate ViperpE,N from EperpN,E using bmag298
eperpe = ipar1d(iepe)*0.01;
ereperpe = ipar1d(ierepe)*0.01;
eperpn = ipar1d(iepn)*0.01;
ereperpn = ipar1d(ierepn)*0.01;
if (ipar1d(iepn) > -32760.)
 vperpe = -ipar1d(iepn)*0.1/bmag298;
else
 vperpe = missing;
end
if (ipar1d(iepe) > -32760.)
 vperpn = ipar1d(iepe)*0.1/bmag298;
else
 vperpn = missing;
end

% Average over j150:nrows for ViperpE,N,par
ne = 0;
nn = 0;
np = 0;
avvpe = 0.;
avvpn = 0.;
avvpp = 0.;
avervpe = 0.;
avervpn = 0.;
avervpp = 0.;
  for j=j150:nrows
%  Delete if Vi error bar > 100 m/s
         if (ipar2d(jervpe,j) > 100. | ipar2d(jvpe,j) < -32760.)
	   ipar2d(jvpe,j) = missing;
	   ipar2d(jervpe,j) = missing;
	 else
	   avvpe = avvpe + ipar2d(jvpe,j);
	   ne = ne + 1;
	   avervpe = avervpe + ipar2d(jervpe,j);
         end
         if (ipar2d(jervpn,j) > 100. | ipar2d(jvpn,j) < -32760.)
	   ipar2d(jvpn,j) = missing;
	   ipar2d(jervpn,j) = missing;
	 else
	   avvpn = avvpn + ipar2d(jvpn,j);
	   nn = nn + 1;
	   avervpn = avervpn + ipar2d(jervpn,j);
         end
         if (ipar2d(jervpp,j) > 100. | ipar2d(jvpp,j) < -32760.)
	   ipar2d(jvpp,j) = missing;
	   ipar2d(jervpp,j) = missing;
	 else
	   avvpp = avvpp + ipar2d(jvpp,j);
	   np = np + 1;
	   avervpp = avervpp + ipar2d(jervpp,j);
         end
  end
if (ne > 0) 
 avvpe = avvpe/ne;
 avervpe = avervpe/ne;
else
 avvpe = missing;
 avervpe = missing;
end
if (nn > 0) 
 avvpn = avvpn/nn;
 avervpn = avervpn/nn;
else
 avvpn = missing;
 avervpn = missing;
end
if (np > 0) 
 avvpp = avvpp/np;
 avervpp = avervpp/np;
else
 avvpp = missing;
 avervpp = missing;
end


fprintf (OUTIW1,' %4d %4d %4d %5.2f %6.2f %6.2f %8.1f %8.1f %8.1f %8.1f %8.1f %8.1f %7.4f %9.2f %9.2f %8.1f %9.2f %9.2f %8.1f\n',ibyrt,ndayn,ibdtt,uth,glat,glon, avvpe,avervpe,avvpn,avervpn,avvpp,avervpp, bmag298,eperpe,ereperpe,vperpn,eperpn,ereperpn,vperpe);


%  Gather data for medians every 20 min at 00,20,40 min
if (imn >= 50 | imn < 10) mm(k20) = 1; end
if (imn >= 10 & imn < 30) mm(k20) = 2; end
if (imn >= 30 & imn < 50) mm(k20) = 3; end
if (mm(k20) == mp(k20))
  nm(k20) = nm(k20) + 1;
else
%  Calculate medians from nm values at mp, then reset nm
%  Consider only for a particular k20 (1 or 2)
	nmk = nm(k20);
clear b20ut;
for n=1:nmk
 b20ut(n) = a20ut(n,k20);
end
	avut = median(b20ut)
	mdayn = a20day(nmk,k20);
	moda = a20moda(nmk,k20);
	if (a20day(1,k20) < a20day(nmk,k20)) 
	  for n=1:nmk-1
	   if (a20day(n,k20) < a20day(nmk,k20)) 
	     b20ut(n) = b20ut(n)-24.;
           end
          end
          avut = median(b20ut);
	 if (avut < 0.) 
	  avut = avut + 24.;
	  mdayn = a20day(1,k20);
	  moda = a20moda(1,k20);	
	 end
        end

for m=1:3
clear var;
clear ervar;
	  nav = 0;
	  for n=1:nmk
	    if (a20vp(n,k20,m) > -32760.) 
	      nav = nav + 1;
	      var(nav) = a20vp(n,k20,m);
	      ervar(nav) = a20ervp(n,k20,m);
	    end
          end
	  if (nav > 0) 
	    avvp20(m) = median(var);
	    avervp20(m) = std(var);
% Save original error bar if std=0 (nav=1)
	    if (nav==1), avervp20(m) = ervar(1); end
	  else
	    avvp20(m) = missing;
	    avervp20(m) = missing;
	  end
          numvp20(m) = nav;
end
for m=1:2
clear var;
clear ervar;
clear varp;
	  nav = 0;
	  for n=1:nmk
	    if (a20ep(n,k20,m) > -327.60) 
	      nav = nav + 1;
	      var(nav) = a20ep(n,k20,m);
	      ervar(nav) = a20erep(n,k20,m);
	      varp(nav) = a20exb(n,k20,m);
	    end
          end
	  if (nav > 0) 
	    avep20(m) = median(var);
	    averep20(m) = std(var);
	    avexb20(m) = median(varp);
	    averexb20(m) = std(varp);
% Save original error bar if std=0 (nav=1)
	    if (nav==1), averep20(m) = ervar(1); end
	  else
	    avep20(m) = missing;
	    averep20(m) = missing;
	    avexb20(m) = missing;
	    averexb20(m) = missing;
	  end
          numep20(m) = nav;
end

%end

fprintf (OUTI20W1,' %4d %4d %4d %5.2f %6.2f %6.2f %2d %8.1f %8.1f %2d %8.1f %8.1f %2d %8.1f %8.1f %7.4f %2d %9.2f %9.2f %8.1f %8.1f %2d %9.2f %9.2f %8.1f %8.1f\n',ibyrt,mdayn,moda,avut,glat,glon, numvp20(1),avvp20(1),avervp20(1),numvp20(2),avvp20(2),avervp20(2),numvp20(3),avvp20(3),avervp20(3), bmag298,numep20(1),avep20(1),averep20(1),avexb20(1),averexb20(1),numep20(2),avep20(2),averep20(2),avexb20(2),averexb20(2));

% Reset nm, clear arrays to find medians again
nm(k20) = 1;


end	% end of calculating av/median etc

	mp(k20) = mm(k20);
	km = nm(k20);
	a20day(km,k20) = ndayn;
	a20moda(km,k20) = ibdtt;
	a20ut(km,k20) = uth;
	a20vp(km,k20,1) = avvpe;
	a20ervp(km,k20,1) = avervpe;
	a20vp(km,k20,2) = avvpn;
	a20ervp(km,k20,2) = avervpn;
	a20vp(km,k20,3) = avvpp;
	a20ervp(km,k20,3) = avervpp;
	a20ep(km,k20,1) = eperpe;
	a20erep(km,k20,1) = ereperpe;
	a20ep(km,k20,2) = eperpn;
	a20erep(km,k20,2) = ereperpn;
	a20exb(km,k20,2) = vperpe;
	a20exb(km,k20,1) = vperpn;
        end

%end	% for ???

end	% of if prologue for header/catalog or data

end	% of while 

nrecrd


return


%cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
function [NDA] = cvt2dn (ITIME);
%          This converts the standard I.S. format time ITIME(4)
%          yyyy, mmdd, hhmm & csec to day number of year NDA.

      NUMD = [0,31,28,31,30,31,30,31,31,30,31,30];
 
      MMDD=ITIME(2);
      MM=fix(MMDD/100);
      DD=fix(MMDD - MM*100);
 
%          Leap year adjustment
      NUMD(3)=28;
if (mod(ITIME(1),4)==0)
NUMD(3)=29;
end
 
      NDA = 0;
for I=1:MM
 NDA = NDA+NUMD(I);
end
      NDA = NDA+DD;
 
if (NDA < 1 | NDA > 366)
' CVT2DN:  Bogus day number %d computed from %d %d %d %d \n',nda,itime;
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
