% foF2 retrievs from ionogram
% with lorentsian type curve
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);

% initial approximation file
fl_appr=sprintf('%s', pth_par, fil_appr);
par_0=load(fl_appr);
sz_par=size(par_0);
L_par_0=sz_par(1);

appr_out_0=zeros(10000,10);
k_a_out=0;
% 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
for i_day=i_day_st:i_day_fn
    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,'\');
    fr_out=zeros(1440,2)-1e11;
    fr_out(:,1)=[1:1440]';
    if key_iono_pl>0
            a_i=figure;
             f_i=subplot(2,1,1);
           set(f_i,'nextplot','replace');
    end
    % cycle over hours
    for i_hr=i_hr_1:i_hr_2
        s_hr=sprintf('%2.2i', i_hr);
        % cycle over minutes
        for i_mnt=i_mnt_1:i_mnt_2
            s_mnt=sprintf('%2.2i', i_mnt);
            s_tlt=sprintf('%s', s_year, '/', s_mon, '/', s_day, ', UT =', s_hr, ':', s_mnt);
            i_mnt_day=i_hr*60+i_mnt;
            if L_par_0>0
            f_par_day=find(par_0(:,1)==i_year & par_0(:,2)==i_doy);
            L_par_day=length(f_par_day);
            
            if L_par_day>0
                par_day=par_0(f_par_day,:);
            
                
                mnt_par=par_day(:,3)*60+par_day(:,4);
                [d_par_min,i_r_par]=min(abs(mnt_par-i_mnt_day));
                par_c=par_day(i_r_par,5:8);
            else
                par_c=par_0(1,5:8);
            end
            else
                % begiining with a random point
            F1=4;
            DF=2.6;
            k=6.5e-3;
            alf=1.37;
            par_c=[F1,DF,k,alf];
            end
            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);
           % frequency and altitude boundaries
           f_fr=find(fr>=fr_int(1) & fr<=fr_int(2));
           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_h_e=length(h_e);
           matr_e=matr(f_alt,f_fr);
       pwr_t=10.^(matr_e/20);
       pwr=pwr_t';
           if key_iono_pl>0
          
           contourf(h_e,fr_e',pwr);
           tlt=title(s_tlt);
        colorbar
           end
        
[cf_lin_1,A_up_dn_1,i_itr_1]=fit_lin_lor(pwr,h_e,fr_e,par_c,h0,h1_int,fr_w_1,pw_b,pw_rat,fr_w_d,cf_pw_dev);
[cf_lin_2,A_up_dn_2,i_itr_2]=fit_lin_lor(pwr,h_e,fr_e,par_c,h0,h2_int,fr_w_2,pw_b,pw_rat,fr_w_d,cf_pw_dev);
if cf_lin_2(2)>cf_lin_1(2) || cf_lin_2(1)<cf_lin_1(1)
    key_new_par=1;
else
    %par_new=appr_2_lor(pwr,h_e,fr_e,appr_dsc,h0,h1_int,h2_int,fr_w_d,cf_lin_1,cf_lin_2);
     [par_new,h_h, H_22]=appr_2_lor_from_2int(pwr,h_e,fr_e,appr_dsc,h0,h1_int,h2_int,fr_w_d,cf_lin_1,cf_lin_2);  
     H_h=(h_h-h0)/H_22;
     F0=par_new(1);
     DF=par_new(2);
     k=par_new(3);
     alf=par_new(4);
     k_out=k/H_22^alf;
     fr_h=F0+DF*(k*H_h.^alf)./(1+k*H_h.^alf);
     fr_h_km=F0+DF*(k_out*h_h.^alf)./(1+k_out*h_h.^alf);
     subplot(f_i);
     set(f_i,'nextplot','add');
     p_apr=plot(h_h, fr_h, 'y');
     set(p_apr,'linewidth', 2);
     p_a2=plot(h_h, fr_h, 'r');
     is_ok=input('Input 1, if you adopt the approximation results   ');
     if is_ok>0
         k_a_out=k_a_out+1;
         foF2=F0+DF;
         k_out=k/H_22^alf;
         appr_out_0(k_a_out,:)=[i_year,i_mon,i_day,i_hr,i_mnt,foF2,F0,DF,k_out,alf];
     end
end
    
    
            end
        end
    end

    
end
end

% final output
if k_a_out>0
    s_mon_1=sprintf('%2.2i', i_mon_1);
    s_mon_2=sprintf('%2.2i', i_mon_2);
    s_day_1=sprintf('%2.2i', i_day_1);
    s_day_2=sprintf('%2.2i', i_day_2);
    appr_out=appr_out_0(1:k_a_out,:);
    fil_appr_out=sprintf('%s', 'sod_', s_year,'_',s_mon_1,'_',s_day_1,'_',s_mon_2,'_',s_day_2,'_foF2_appr.txt');
    fil_appr_lst=sprintf('%s', 'sod_', s_year,'_',s_mon_1,'_',s_day_1,'_',s_mon_2,'_',s_day_2,'_foF2_appr.lst');
    fl_appr_out=sprintf('%s',pth_out_b,fil_appr_out);
    fl_appr_lst=sprintf('%s',pth_out_b,fil_appr_lst);
    save(fl_appr_out,'appr_out','-ascii');
    i_f_lst=fopen(fl_appr_lst,'w');
    fprintf(i_f_lst, '%s\n', 'Columns description');
    fprintf(i_f_lst, '%s\n', '___________________');
    fprintf(i_f_lst, '%s\n', '1. Year');
    fprintf(i_f_lst, '%s\n', '2. Month');
    fprintf(i_f_lst, '%s\n', '3. Day');
    fprintf(i_f_lst, '%s\n', '4. Hour');
    fprintf(i_f_lst, '%s\n', '5. Minute');
    fprintf(i_f_lst, '%s\n', '6. foF2, MHz');
    fprintf(i_f_lst, '%s\n', '___________________');
    fprintf(i_f_lst, '%s\n', 'F(h)=F0+DF*(k*(h-h0)^alf)/(1+k*(h-h0)^alf)');
    fprintf(i_f_lst, '%s\n', '___________________');
    fprintf(i_f_lst, '%s\n', '7. F0');
    fprintf(i_f_lst, '%s\n', '8. DF');
    fprintf(i_f_lst, '%s\n', '9. k');
    fprintf(i_f_lst, '%s\n', '10. alf');
    fclose(i_f_lst);
    
    
end