% selects critical frequency from a difital ionogram

clear all;
close all;
pth_par='par\';
fil_mask='*.par';
mask_full=sprintf('%s', pth_par, fil_mask);
[fil_par, pathname] = uigetfile(mask_full, 'Input parameter file');
[vrb,blk]=rd_par_f(pth_par,fil_par);
s_v=size(vrb);
l_v=s_v(1,1); % length
% assignment of variables
for i_v=1:l_v
    assignin('base', vrb{i_v,2},vrb{i_v,3}); % Asignment of variables
end
% assignment of cell arrays
s_b=size(blk);
l_b=s_b(1,1);
for i_b=1:l_b
    s_in=size(blk{i_b,3});
    w_in=s_in(1,2);
    assignin('base', blk{i_b,2},blk{i_b,3});
end
s_year=sprintf('%4.4i', i_year);
% maximal altitude
h_max_layr=alt_int(2);

% cycle over days
for i_mon=i_mon_1:i_mon_2
    s_mon=sprintf('%2.2i', i_mon);
    if i_mon==i_mon_1
        i_day_st=i_day_1;    
    else
        i_day_st=1;
    end
    if i_mon==i_mon_2
        i_day_fn=i_day_2;    
    else
        i_day_fn=mon_lngth(i_year,i_mon);
    end
k_mnt=0;
pth_out=sprintf('%s',pth_out_b,s_year,'\');
pth_fig=sprintf('%s',pth_fig_b,'mlab_fig\',s_year,'\');
pth_jpg=sprintf('%s',pth_fig_b,'jpeg\',s_year,'\');
for i_day=i_day_st:i_day_fn
    i_doy=day_of_year(i_year,i_mon,i_day);
    
    s_day=sprintf('%2.2i', i_day);
    i_doy=day_of_year(i_year,i_mon,i_day);
    s_doy=sprintf('%3.3i', i_doy);
    pth_in=sprintf('%s',pth_in_b,sens,s_year,s_mon,s_day,'\');
    L_mnt_day=(i_hr_2-i_hr_1+1)*60;
    fr_out=zeros(L_mnt_day,2)-1e11;
    fr_out(:,1)=[i_hr_1*3600:60:(i_hr_2+1)*3600-60]';
    k_mnt_day=0;
    fil_out=sprintf('%s', 'sgo_',s_year,'_',s_doy,'_',fil_out_a,'.txt');
       fil_fig=sprintf('%s', 'sgo_',s_year,'_',s_doy,'_',fil_out_a,'.fig');
    fil_jpg=sprintf('%s', 'sgo_',s_year,'_',s_doy,'_',fil_out_a,'.jpeg');
    fl_out=sprintf('%s',pth_out,fil_out);
    if key_iono_pl>0
            a_i=figure;
             
    end
    % cycle over hours
    for i_hr=i_hr_1:i_hr_2
        s_hr=sprintf('%2.2i', i_hr);
        % cycle over minutes
        if i_hr==i_hr_1
            i_mnt_st=i_mnt_1;
        else
            i_mnt_st=0;
        end
        if i_hr==i_hr_2
            i_mnt_fn=i_mnt_2;
        else
            i_mnt_fn=59;
        end
        for i_mnt=i_mnt_st:i_mnt_fn
            dat_tim=[i_year,i_doy,i_hr,i_mnt];
k_mnt=k_mnt+1;
k_mnt_day=k_mnt_day+1;
h_fr_c=[];
            i_mnt_day=i_hr*60+i_mnt;
            i_mnt_c=i_mnt_day+1;
            s_mnt=sprintf('%2.2i', i_mnt);
            fil_in=sprintf('%s', sens,s_year,s_mon,s_day,'_',s_hr,s_mnt);
            fil_in_cntrl=sprintf('%s', sens,s_year,s_mon,s_day,'_',s_hr,s_mnt,'.mat');
            f_in=sprintf('%s', pth_in,fil_in);
            fl_in_cntrl=sprintf('%s', pth_in,fil_in_cntrl);
            i_f_in=fopen(fl_in_cntrl);
            if i_f_in>0
                fclose(i_f_in);
            f1=load(f_in);
            ionogr=f1.ionogr;
            h_lim=ionogr.ylim;
            fr_lim=ionogr.xlim;
            matr=ionogr.data;
            s_matr=size(matr);
            L_fr=s_matr(2);
            L_h=s_matr(1);
           fr=linspace(fr_lim(1),fr_lim(2),L_fr);
           h=linspace(h_lim(1),h_lim(2),L_h);
           d_fr=(fr_lim(2)-fr_lim(1))/(L_fr-1);
           % frequency and altitude boundaries
           f_fr=find(fr>=fr_int(1) & fr<=fr_int(2));
           f_fr_0=find(fr>=fr_int(1)-fr_hw & fr<=fr_int(2)+fr_hw);
           L_fr=length(f_fr);
           f_alt=find(h>=alt_int(1) & h<=alt_int(2));
           fr_e=fr(f_fr);
           h_e=h(f_alt);
           L_alt=length(f_alt);
           fr_ee=ones(L_alt,1)*fr_e;
           L_h_e=length(h_e);
           matr_e0=matr(f_alt,f_fr_0);
           fr_fr=fr(f_fr);
           matr_e=zeros(L_h_e,L_fr);
           for i_fr=1:L_fr
               f_fr_c=find(fr_fr>=fr_fr(i_fr)-fr_hw & fr_fr<fr_fr(i_fr)+fr_hw);
               L_fr_c=length(f_fr_c);
               if L_fr_c>0
               matr_c=matr_e0(:,f_fr_c);
               matr_e(:,i_fr)=(mean(matr_c'))';
               end
           end
           %matr_e=matr(f_alt,f_fr);
           if key_iono_lin>0
       pwr=10.^(matr_e/20);
           else
               pwr=matr_e;
           end
           if key_iono_pl>0
          
          % contourf(fr_e,h_e',pwr);
          f_i=subplot(1,1,1);
          set(f_i,'nextplot','replace'); 
contourf(h_e,fr_e',pwr');
        colorbar
        s_tltt=sprintf('%s', s_year,', DOY=', s_doy, ', UT=', s_hr,':',s_mnt);
    tltt=title(s_tltt);
    set(f_i,'nextplot','add'); 
        end
           
% 1. According to maximal/mean power in an altitude band
if key_alg==1
    f_hb=find(h_e>=alt_b(1) & h_e<=alt_b(2));
    pwr_b=pwr(f_hb,:);
        % average power in altitude band
    if key_av==1
        pwr_ref=mean(pwr_b');
    elseif key_av==3
        pwr_ref=max(pwr_b');
    end
    key_new=1;
    for i_fr=fr_hw+2:L_fr-fr_hw
        pwr_left=pwr_ref(i_fr-fr_hw-1:i_fr-1);
        pwr_right=pwr_ref(i_fr:i_fr+fr_hw);
    %pause
       if key_av==1
        pwr_lf=mean(pwr_left);
        pwr_rt=mean(pwr_right);
       % pause
    elseif key_av==3
        pwr_lf=max(pwr_left);
        pwr_rt=max(pwr_right);
       end 
    if (pwr_lf>pwr_thrsh && pwr_rt<pwr_thrsh*cf_thrsh);
        
        fr_c=fr_e(i_fr);
        key_new=0;
    %    pause
        fr_out(i_mnt_c,2)=fr_e(i_fr);
    end
    end
elseif  key_alg==2
    i_c=0;
 fr_pwr=[]; 
 fr_prev=0;
 dfr_prev=0;
    for i_h=h_hw+1:L_h_e-h_hw
        i_c=i_c+1;
        pwr_c=pwr(i_h-h_hw:i_h+h_hw,:);
        pw_c=mean(pwr_c);
       % if key_iono_pl>0
        %     f_i2=subplot(2,1,2);
          
         %  plot(fr_e,pw_c);
        %end
        [fr_pw]=find_lmax_iono(pw_c,fr_e,fr_int,pwr_mm,num_bands_in);
       s_frmx=size(fr_pw);
        L_frmx=s_frmx(1);
        fr_mx(i_c,:)=zeros(1,num_bands_in)-1e11;
        fr_mx(i_c,1:L_frmx)=fr_pw(:,1)';
        pw_mx(i_c,:)=zeros(1,num_bands_in)-1e11;
        pw_mx(i_c,1:L_frmx)=fr_pw(:,2)';
        h_mx(i_c,1)=h_e(i_h);
        n_mx(i_c,:)=zeros(1,num_bands_in)-1e11;
        n_mx(i_c,1:L_frmx)=fr_pw(:,3)';
    end
  % [h_gr,fr_gr,hfp_gr]= pl_pw_iono(h_mx,fr_mx,pw_mx,dh_pl_dsc,dfr_pl_dsc);
   %[h_b,fr_b]=fr_bound_iono(h_gr,fr_gr,hfp_gr,band_dsc);
   fil_prm_mult=sprintf('%s', fil_prm_mult_b, '_', fil_prm_mult_a, '.', extn_prm_mult);
   fil_prm_out=sprintf('%s', fil_prm_mult_b, '_', fil_prm_mult_a, '.', extn_prm_out);
   %[par_out,rho_out]=appr_lor_iono_line(h_mx,fr_mx,pw_mx,h_appr_dsc,pth_par,fil_prm_mult,fil_prm_out,lin_dsc,ar_dsc,dlt_par,eps_par,cf_pw_rat,dat_tim,bcgr_dsc); 
   [par_out,rho_out]=appr_lor_iono_line(h_e',fr_ee,pwr,h_appr_dsc,pth_par,fil_prm_mult,fil_prm_out,lin_dsc,ar_dsc,dlt_par,eps_par,cf_pw_rat,dat_tim,bcgr_dsc); 
  
   h1=h_appr_dsc(1);
    h2=h_appr_dsc(2);
    h0=h_appr_dsc(3);

f_h_appr=find(h_mx>=h1 & h_mx<=h2);
L_h_appr=length(f_h_appr);
if L_h_appr>0
    ha=h_mx(f_h_appr);
   fr_c=par_out(1);
        df_c=par_out(2);
        k_c=par_out(3);
        alf_c=par_out(4);
        ya=fr_c+df_c*(k_c*(ha-h0).^alf_c./(k_c*(ha-h0).^alf_c+1));
        par_out_min=min(par_out);
        if par_out_min>0  && key_iono_pl>0
            ppapl=plot(ha,ya,'w');
            set(ppapl, 'linewidth', 2);
        end
end
   
  if key_oth>0  % drafts of previous versions   
    fl_prm_mult=sprintf('%s', pth_par, fil_prm_mult);
    i_f_prm_mult=fopen(fl_prm_mult);
    if i_f_prm_mult>0
        fclose(i_f_prm_mult);
   prm_mult=load(fl_prm_mult);
 [prm_out,dsc_out]=fr_crit_mult_lor(n_mx,fr_mx,h_mx,pw_mx,h_int,num_bands_in,num_bands_out,prm_mult,dp_mult,eps_mult,cf_alf_k,lin_dsc,ar_dsc);   
 [pw_max,ih_pw_max]=max(pw_mx);
 [pw_max_gl,ifr_pw_mx]=max(pw_max);
 ih_pw_max_gl=ih_pw_max(ifr_pw_mx);
 ifr_max_gl=n_mx(ih_pw_max_gl,ifr_pw_mx);
 fr_max_gl=fr_e(ifr_max_gl);
 ih_max_gl=ih_pw_max_gl-h_hw;
 if ih_pw_max_gl<=1  
     ih_pw_max_gl=2;
 end
 %h_max_gl=h_e(ih_max_gl);
 h_max_gl=h_e(ih_pw_max_gl);
 h_max_gl_0=h_e(ih_pw_max_gl-1);
 pwr_1=pwr(ih_pw_max_gl,:);
 pwr_0=pwr(ih_pw_max_gl-1,:);
 np_fr_gap=floor(w_fr_gap/d_fr);
 [i_fr_cr_1,fr_cr_1,df_dh_10,d2f_dh_10]=frq_crit(ih_pw_max,pwr_1,fr_e,h_e,ifr_max_gl,np_fr_gap,key_iono_lin,pw_rat,key_av_bcgr,0,0,0,0,d_frc_max);
 [i_fr_cr_0,fr_cr_0,df_fh_00,d2f_dh_00]=frq_crit(ih_pw_max-1,pwr_0,fr_e,h_e,ifr_max_gl,np_fr_gap,key_iono_lin,pw_rat,key_av_bcgr,0,0,0,0,d_frc_max);
 df_dh_1=(i_fr_cr_1-i_fr_cr_0)/(h_max_gl-h_max_gl_0);
 k_h=1;
 h_fr_c(k_h,:)=[ih_pw_max_gl,h_max_gl,i_fr_cr_1,fr_cr_1,df_dh_1,d2f_dh_10];
 while h_max_gl>h_max_layr
     % to add if the initial altitude is too big
 end
 ih_prev=ih_max_gl;
 ih_pr2=ih_max_gl-1;
 h_prev=h_max_gl;
 fr_max_prev=fr_max_gl;
 ifr_max_prev=ifr_max_gl;
 fr_cr_prev=fr_cr_1;
 i_fr_cr_prev=i_fr_cr_1;
 df_dh_prev=df_dh_1;
 if ih_max_gl<1
     ih_max_gl=1;
 end
 %for ih=41:45
 for ih=ih_max_gl+1:L_h_e-h_hw
     ic=ih-h_hw;
     a_fr_diff_c=abs(fr_mx(i_c,:)-fr_max_prev);
     pwr_c=pwr(ih,:);
   %  [adfr_min,i_adfr_min]=min(a_fr_diff_c);
     %if adfr_min<=w_fr_gap
         %ifr_max_c=n_mx(i_c,i_adfr_min);
         %fr_max_c=fr_mx(i_c,i_adfr_min);
         [pw_max_c,ifr_max_c]=max(pwr_c);
         fr_max_c=fr_e(ifr_max_c);
         [i_fr_cr_c,fr_cr_c,df_dh_c,d2f_dh_c]=frq_crit(ih,pwr_c,fr_e,h_e,ifr_max_c,np_fr_gap,key_iono_lin,pw_rat,key_av_bcgr,i_fr_cr_prev,ih_prev,ih_pr2,df_dh_prev,d_frc_max);
         if abs(fr_cr_c-fr_cr_prev)<=w_fr_gap
             k_h=k_h+1;
             h_fr_c(k_h,:)=[ih,h_e(ih),i_fr_cr_c,fr_cr_c,df_dh_c,d2f_dh_c];
             ih_pr2=ih_prev;
             ih_prev=ih;
             
 h_prev=h_e(ih);
 fr_max_prev=fr_max_c;
 ifr_max_prev=ifr_max_c;
 fr_cr_prev=fr_cr_c;
 i_fr_cr_prev=i_fr_cr_c;
         else
             % to add if critical frequencies in adjacent layers  differ
         % essentially 
         end
         %else
         % to add if maximal frequencies in adjacent layers  differ
         % essentially 
  %   end
 end
     
 
 
        L_frmx=length(fr_mx); 
        
        if L_frmx>0
        
        if i_c>=2
            dfr_prev=fr_mx-fr_prev;
        end
        fr_prev=fr_mx;
        %fr_mx=fr_e(i_fr);
      %  h_pwr_0(i_c,:)=[fr_mx,h_mx,pw_mx];
        
      
        end
    end
    %f_ok_pwr=find(h_pwr_0(:,3)>0);
    %h_pwr=h_pwr_0(f_ok_pwr,:);
%plot(h_pwr(:,1),h_pwr(:,3));
%set(f_g, 'xlim', [0 3],'xgrid', 'on');
%df_c=-diff(h_pwr(:,3))./diff(h_pwr(:,1));
%[m_dfc,im_dfc]=max(df_c);
%fr_out(i_mnt_c,2)=fr_e(im_dfc);
%fr_out(i_mnt_c,2)=h_pwr(im_dfc,1);
%pause
        
        

        
        

if k_mnt>1 && alf_a>0  && k_a>0
alf_prev=alf_a;
k_prev=k_a;
fr_int_prev=fr_int_a;
else
alf_prev=1.5;
k_prev=.003;
fr_int_prev=[5.7,8.5];
end
if key_iono_pl>0
subplot(f_i);
set(f_i,'nextplot','add');
pa=plot(h_fr_c(:,2),h_fr_c(:,4),'r--');
set(pa,'linewidth',3);
end
[alf_a,k_a,fr_int_a,da_alf,da_k,dqa_ea]=appr_lor(h_fr_c,h_int,alf_prev,k_prev,fr_int_prev,k_mnt,dp_mult,eps_mult,num_steps_max,cf_alf_k,key_both);
if alf_a>0 && k_a>0
f_hf_pl=find(h_fr_c(:,2)>h_int(1) & h_fr_c(:,2)<h_int(2));
hf_pl=h_fr_c(f_hf_pl,[2,4,3]);
fr_appr=fr_int_a(1)+fr_int_a(2)*k_a*(hf_pl(:,1)-h_int(1)).^alf_a./(k_a*(hf_pl(:,1)-h_int(1)).^alf_a+1);
L_frapr=length(fr_appr);
if L_frapr>0
hf_pl(:,3)=fr_appr;
if key_iono_pl>0
p_apr=plot(hf_pl(:,1),hf_pl(:,3),'m');
        set(p_apr,'linewidth',2);
end
        fr_crt_c=fr_appr(L_frapr);
        
    fr_out(k_mnt_day,2)=fr_crt_c;
end
end
        end
            end
        end
        end
    end
    fl_prm_out=sprintf('%s', pth_par,fil_prm_out);
    i_f_prm=fopen(fl_prm_out);
    if i_f_prm>0
        fclose(i_f_prm);
        prm_pl=load(fl_prm_out);
        aap=figure;
        t_prm=prm_pl(:,3)+prm_pl(:,4)/60;
        ts_prm=prm_pl(:,3)*3600+prm_pl(:,4)*60;
        h_ap=h_int(2)-h_appr_dsc(3);
        fr_appr=prm_pl(:,5)+prm_pl(:,6).*(prm_pl(:,7).* h_ap.^prm_pl(:,8)./(1+ prm_pl(:,7).*h_ap.^prm_pl(:,8)));
        f0f2_out=[ts_prm,fr_appr];
        f_ap=subplot(1,1,1);
        p_ap=plot(t_prm,fr_appr);
        set(f_ap,'fontname', 'Times New Roman', 'fontsize', 14,'xgrid', 'on','ygrid','on');
        set(p_ap,'color','b','linewidth',1.5);
        y_tlt=ylabel('Fof2, MHz');
        x_tlt=xlabel('UT, hours');
        tlt_txt=sprintf('%s', s_year,', DOY=', s_doy);
        tlt_fig=title(tlt_txt);
        fl_fig=sprintf('%s', pth_fig,fil_fig);
       fl_jpg=sprintf('%s', pth_jpg,fil_jpg);
       saveas(aap,fl_fig,'fig');
       saveas(aap,fl_jpg,'jpg');
       save(fl_out, 'f0f2_out', '-ascii');
    end
   %save(fl_out,'fr_out','-ascii');     
   %f_ok_d=find(fr_out(:,2)>0);
   %L_ok_d=length(f_ok_d);
   %if L_ok_d>0
   %    fr_ok=fr_out(f_ok_d,:);
   %    aad=figure;
   %    fdd=subplot(1,1,1);
   %    plot(fr_ok(:,1)/3600,fr_ok(:,2));
   %    tlt_txt=sprintf('%s', s_year,', DOY=', s_doy);
   %    tlt_f=title(tlt_txt);
   %    fl_fig=sprintf('%s', pth_fig,fil_fig);
   %    fl_jpg=sprintf('%s', pth_jpg,fil_jpg);
   %    saveas(aad,fl_fig,'fig');
   %    saveas(aad,fl_jpg,'jpg');
   
   
    %end
    
end
end



