
function out = rd_son
%%
%% This program is written for reading Sonde ISR CEDAR DB .asc data
%% BA Emery, April 2011
clear all; close all
warning off MATLAB:divideByZero
 
%  kindat=5501 long-pulse ~4.5min 3-position (Z, NE+NW 75elev) 98-1063/95-1053km
%    Pick every Z position (cannot average w 75elev since at diff alts) from 98-1063km
%  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;
% ifc75 = 1, save also the 75 deg elevation data (2 pos to NE and NW w 1 more ht than htz)
	ifc75 = 0

% Open 2 files - one for 1D TEC, hmax, Nmax, and FoF2, and the other for 2D ht profiles
%  Make 2 extra files of medians every 20 min (4 if save 75 deg elev also)

% Long Pulse (LP) with Z, NE, NW, ..., Z, NE (no NW at end); ifc75=1, use NE time
% Set previous minute to same as first time read
	mp(1) = 3;
	mp(2) = 3;
%  datfil = sprintf ('son050514g.003_5501_98-1063km');
% iw1 = sprintf ('son_lp050514_15_nmax');
% iw2 = sprintf ('son_lp050514_15_profiles');
% i20w1 = sprintf ('son_lp050514_15_20m_nmax');
% i20w2 = sprintf ('son_lp050514_15_20m_profiles');
	mp(1) = 2;
	mp(2) = 2;
  datfil = sprintf ('son061214g.007_22_1215_22UT_3pos_5501');
  iw1 = sprintf ('son_lp061214_15_nmax');
  iw2 = sprintf ('son_lp061214_15_profiles');
  i20w1 = sprintf ('son_lp061214_15_20m_nmax');
  i20w2 = sprintf ('son_lp061214_15_20m_profiles');
  DAT = fopen (datfil,'r');

	nm(1) = 0;
	nm(2) = 0;
	glat = 66.99;
	glon = -50.95;
%  Ne,Ti,Te,LOS and errs at 11,15,13,17 12,16,14,18
	jne = 11;
	jerne = 12;
	jti = 15;
	jerti = 16;
	jte = 13;
	jerte = 14;
	jvi = 17;
	jervi = 18;
% 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 99 and 120 km for long-pulse (or 95 and 116 km for 75 deg elev)
	j150 = 2;
	j150 = 3;


  OUTIW1 = fopen (iw1,'w');
 fprintf (OUTIW1,' year dayn moda   uth   glat   glon   hmax foF2MHz  Nmax(m-3)    TEC(m-2) \n');

  OUTIW2 = fopen (iw2,'w');
 fprintf (OUTIW2,' year dayn  uth  alt(km) Vi m/s  errVi  Ti(K)  errTi  Te(K)  errTe lg10Ne(m-3)   errNe \n');

  OUTI20W1 = fopen (i20w1,'w');
 fprintf (OUTI20W1,' year dayn moda   uth   glat   glon   hmax foF2MHz  Nmax(m-3)    TEC(m-2) \n');

  OUTI20W2 = fopen (i20w2,'w');
 fprintf (OUTI20W2,' year dayn   uth alt(km) #Vi medViz(m/s) stdVi #Ti medTi(K)    stdTi #Te medTe(K)    stdTe #Ne   medNe(m-3)   stdNe \n');
if (ifc75==1)
  k20w1 = sprintf ('son_lp75_050514_15_20m_nmax');
  OUTK20W1 = fopen (k20w1,'w');
 fprintf (OUTK20W1,' year dayn moda   uth   glat   glon   hmax foF2MHz  Nmax(m-3)    TEC(m-2) \n');
  k20w2 = sprintf ('son_lp75_050514_15_20m_profiles');
  OUTK20W2 = fopen (k20w2,'w');
 fprintf (OUTK20W2,' year dayn   uth alt(km) #Ti medTi(K)    stdTi #Te medTe(K)    stdTe #Ne   medNe(m-3)   stdNe \n');
end

        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

%  Find if ipar1d(2) > 8990 (elev*100)

elev = ipar1d(2)*0.01;

% Consider hmax only below 600 km; LP htz(13) = 560km, htjm(14)=563km for 75 deg elev
% 22 Apr 2011:  Anja Stromme suggests cutting all off at 500 km (Z#12 497km, 75d#13 502 km)
j600 = 13;
k20 = 2;
if (elev > 89.)
 ifrd = 1;
 j600 = 12;
 k20 = 1;
 if (nrecrd < 4) 
  nhtz = nrows;
 end
elseif (ifc75==1)
 ifrd = 2;
else
 ifrd = 0;
end
j600m1 = j600 -1;

if (ifrd > 0)

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

clear htmj;
clear corrne;
clear ercorrne;

% 
  for j=1:nrows
%  Set htjm altitude and from first Zenith, set htz altitude (75 elev 1 more ht)
    if (nrecrd < 4 & ifrd == 1)
       htz(j) = ipar2d(1,j)+ipar2d(2,j)*0.0001;
    end
	htjm(j) = ipar2d(1,j)+ipar2d(2,j)*0.0001;
%  Delete if Ne error bar > Ne
         if (ipar2d(jne+1,j) > ipar2d(jne,j))
	   ipar2d(jne,j) = missing;
	   ipar2d(jne+1,j) = missing;
         end
%	rnelog10(j) = ipar2d(jne,j);
%  Delete if Te, Ti error bars > 200K
         if (ipar2d(jerti,j) > 200.)
	   ipar2d(jerti,j) = missing;
	   ipar2d(jti,j) = missing;
         end
         if (ipar2d(jerte,j) > 200.)
	   ipar2d(jerte,j) = missing;
	   ipar2d(jte,j) = missing;
         end
%  Delete if Ti < 500K when ht > 300km
         if (ipar2d(jti,j) < 500. & htjm(j) > 300.)
	   ipar2d(jerti,j) = missing;
	   ipar2d(jti,j) = missing;
	 end
%  Delete if Vi LOS error bars > 200 m/s
         if (ipar2d(jervi,j) > 200.)
	   ipar2d(jervi,j) = missing;
	   ipar2d(jvi,j) = missing;
         end
  end

%  Could convert to htz (nhtz) from htjm (nhtz+1) for the 75 deg elev files
%    but prefer to write another code to combine the 20-min averages

  for j=1:nrows
%  Set Ne in m-3
	 corrne(j) = missing;
	 ercorrne(j) = missing;
 	 if (ipar2d(jne,j) > 0.)
	  corrne(j) = 10.^(ipar2d(jne,j)*0.001);
	  ercorrne(j) = 10.^(ipar2d(jne+1,j)*0.001);
	 end
  end
% ipar2d(jne,:)
%  Set top rnelog10 if missing
%  if (rnelog10(nrows) < 0.) 
%   for j=nrows-1:-1:j100p
%	   if (rnelog10(j) > 0. & rnelog10(j-1) > 0. & rnelog10(j-1) > rnelog10(j)) 
%	      diffne = (rnelog10(j-1)-rnelog10(j))/(htjm(j)-htjm(j-1));
%              for j1=j+1:nrows
%	        if (rnelog10(j1) < 0.) 
%                 rnelog10(j1) = rnelog10(j) - diffne*(htjm(j1)-htjm(j));
%	        end
%              end
%	   end
%   end
%  end
% rnelog10(nrows)
%  Set "bottom" rnelog10 if missing
%	if (rnelog10(j100) < 0.)
%	  for j=j100+1:nrows
%	    if (rnelog10(j) > 0.), break, end
%         end
%          j1 = j;
%	  for j=j1:nrows
%	    if (rnelog10(j) > 0.), break, end
%          end
%          j2 = j;
%	  diffne = (rnelog10(j2)-rnelog10(j1))/(htjm(j2)-htjm(j1));
%          rnelog10(j100) = rnelog10(j1) - diffne*(htjm(j1)-htjm(j100));
%    	end
% rnelog10(j100)
%  Fill in missing rnelog10 values as linear interpolations
%	for j=j100+1:nrows-1
%	  if (rnelog10(j) < 0.)
%	    for j1=j,nrows-1
%	      if (rnelog10(j1) > 0.), break, end
%           end
%            j2 = j1;
%	    diffne = (rnelog10(j1)-rnelog10(j-1))/(htjm(j1)-htjm(j-1));
%	    rnelog10(j) = rnelog10(j-1) + diffne*(htjm(j)-htjm(j-1));
%	  end
%       end

%  Save 1-min data only for Zenith profiles
if (ifrd==1)

%  Calculate hmF2, NmF2, and TEC (1 TECU=10^16elec/m2)
	xnemax = 0.;
	for j=j150:j600m1
%  Find max Ne from 154.6-536.8 km
	  if (corrne(j) > xnemax)
	    xnemax = corrne(j);
	    jx = j;
	  end
        end
	h(2) = htjm(jx);
	p(2) = corrne(jx);
	hmf2 = htjm(jx);
	xnmax = corrne(jx);
%  Check to see if max Ne next to missing data or not
	for j=jx-1:-1:j150
	  if (corrne(j) > 0.), break, end
        end
       if (j > j150)
	h(1) = htjm(j);
	p(1) = corrne(j);
	for j=jx+1:j600
	  if (corrne(j) > 0.), break, end
        end
       if (j < j600)
	h(3) = htjm(j);
	p(3) = corrne(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 or from j150-j600
	tec = 0.;
	htjm1 = htjm(max(1,j100m1));
%       for j=j100p:nrows
        for j=j150:j600
	  if (corrne(j) > 0.) 
	    tec = tec + corrne(j)*(htjm(j)-htjm1)*1.e+3;
	    htjm1 = htjm(j);
	  end
        end

fprintf (OUTIW1,' %4d %4d %4d %5.2f %6.2f %6.2f %7.2f %5.2f %11.4e %11.4e \n',ibyrt, ndayn,ibdtt,uth,glat,glon,hmf2,fof2,xnmax,tec);

% Save only for j150-j600 
%	for j=j100:nrows
	for j=j150:j600
fprintf (OUTIW2,' %4d %4d %5.2f %7.2f %6d %6d %6d %6d %6d %6d %7.3f %7.3f \n',ibyrt, ndayn,uth,htjm(j),ipar2d(jvi,j),ipar2d(jervi,j),ipar2d(jti,j),ipar2d(jerti,j),ipar2d(jte,j),ipar2d(jerte,j),ipar2d(jne,j)*0.001,ipar2d(jerne,j)*0.001);  
       end 

end	% ifrd==1 (Zenith data)

%  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 j=1:nrows
clear var;
clear ervar;
	  nav = 0;
	  for n=1:nmk
	    if (a20ne(j,n,k20) > 0.) 
	      nav = nav + 1;
	      var(nav) = a20ne(j,n,k20);
	      ervar(nav) = a20erne(j,n,k20);
	    end
          end
	  if (nav > 0) 
	    avne(j) = median(var);
	    averne(j) = std(var);
% Save original error bar if std=0 (nav=1)
	    if (nav==1), averne(j) = ervar(1); end
	  else
	    avne(j) = missing;
	    averne(j) = missing;
	  end
          numne = nav;
clear var;
clear ervar;
	  nav = 0;
	  for n=1:nmk
	    if (a20te(j,n,k20) > 0.) 
	      nav = nav + 1;
	      var(nav) = a20te(j,n,k20);
	      ervar(nav) = a20erte(j,n,k20);
	    end
          end
	  if (nav > 0) 
	    avte(j) = median(var);
	    averte(j) = std(var);
% Save original error bar if std=0 (nav=1)
	    if (nav==1), averte(j) = ervar(1); end
	  else
	    avte(j) = missing;
	    averte(j) = missing;
	  end
          numte = nav;
clear var;
clear ervar;
	  nav = 0;
	  for n=1:nmk
	    if (a20ti(j,n,k20) > 0.) 
	      nav = nav + 1;
	      var(nav) = a20ti(j,n,k20);
	      ervar(nav) = a20erti(j,n,k20);
            end
          end
	  if (nav > 0) 
	    avti(j) = median(var);
	    averti(j) = std(var);
% Save original error bar if std=0 (nav=1)
	    if (nav==1), averti(j) = ervar(1); end
	  else
	    avti(j) = missing;
	    averti(j) = missing;
	  end
          numti = nav;
clear var;
clear ervar;
	  nav = 0;
	  for n=1:nmk
	    if (a20vi(j,n,k20) > 0.) 
	      nav = nav + 1;
	      var(nav) = a20vi(j,n,k20);
	      ervar(nav) = a20ervi(j,n,k20);
            end
          end
	  if (nav > 0) 
	    avvi(j) = median(var);
	    avervi(j) = std(var);
% Save original error bar if std=0 (nav=1)
	    if (nav==1), avervi(j) = ervar(1); end
	  else
	    avvi(j) = missing;
	    avervi(j) = missing;
	  end
          numvi = nav;
% Print out j150-j600
%	  if (j >= j100)
	if (j >= j150 & j <= j600)
if (k20==1)
fprintf (OUTI20W2,' %4d %4d %5.2f %7.2f %3d %8.1f %8.1f %3d %8.1f %8.1f %3d %8.1f %8.1f %3d %11.4e %11.4e \n',ibyrt,mdayn,avut,htjm(j),numvi,avvi(j),avervi(j),numti,avti(j),averti(j),numte,avte(j),averte(j),numne,avne(j),averne(j) ); 
else
fprintf (OUTK20W2,' %4d %4d %5.2f %7.2f %3d %8.1f %8.1f %3d %8.1f %8.1f %3d %8.1f %8.1f %3d %11.4e %11.4e \n',ibyrt,mdayn,avut,htjm(j),numvi,avvi(j),avervi(j),numti,avti(j),averti(j),numte,avte(j),averte(j),numne,avne(j),averne(j) ); 
end
          end   
end
%  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:nrows
        for j=j150:j600
	  if (avne(j) > 0.) 
	    tec = tec + avne(j)*(htjm(j)-htjm1)*1.e+3;
	    htjm1 = htjm(j);
	  end
        end
if (k20==1)
fprintf (OUTI20W1,' %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);
else
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

% 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;
	for j=1:nrows
	  a20ne(j,km,k20) = corrne(j);
	  a20erne(j,km,k20) = ercorrne(j);
	  a20te(j,km,k20) = ipar2d(jte,j);
	  a20erte(j,km,k20) = ipar2d(jerte,j);
	  a20ti(j,km,k20) = ipar2d(jti,j);
	  a20erti(j,km,k20) = ipar2d(jerti,j);
	  a20vi(j,km,k20) = ipar2d(jvi,j);
	  a20ervi(j,km,k20) = ipar2d(jervi,j);
        end

end	% for ???

end	% end of if elev > 89 (Zenith) or ifrd>0

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
