close all;
clear all;

r_s=6.957e10;
r1au=1.496e13;

gauss2nt=1e5;

data_dir = '../data/';
results_dir = '../results/';
%
% Diachronic observed CHM map from database:
% The first file is the interpolated map (floating point)
% The second file is the ceiling of the first (binary furthest extent)
% which is used to identify the hole boundaries to integrate the interpolated map within.
%
chm_synth_file_name        = [data_dir 'CHM/' 'CR2098_D_VSM_PFIT_MHD_CHMSYNTH.hdf'];
chm_synth_binary_file_name = [data_dir 'CHM/' 'CR2098_D_VSM_PFIT_MHD_CHMSYNTH_BINARY.hdf'];

map_root_vec={ ...
'CR2098_D_VSM_PFIT' ...
};

output_file_name='ofp_tableX_mhd_synthchm_flux_results';
[~,~]=system(['rm ' results_dir output_file_name '.txt']);

diary([results_dir output_file_name '.txt']);

disp([sprintf('%30s ','MAP NAME'),...
      sprintf('%12s ','Area'),...
      sprintf('%12s ','UnsignedFlux'),...
      sprintf('%12s ','OpenFlux'),...
      sprintf('%12s ','FluxBalErr%'),...    
      sprintf('%12s ','Br1AU'),...          
     ]);  

idx_start=1;
idx_end=length(map_root_vec);

% Load observed CHM:
[CHM,tvec,pvec] = rdhdf_2d(chm_synth_file_name);
CHM = double(CHM);

% Load observed CHM binary extent:
[CHM_BINARY,~,~] = rdhdf_2d(chm_synth_binary_file_name);
CHM_BINARY = double(CHM_BINARY);

% Get invidual holes' idx from binary map
[CHMS, num_CHMS] = bwlabel(CHM_BINARY, 4);

% Get 2D scale information and sint:
tvec = double(tvec);
pvec = double(pvec);
[~,T]=meshgrid(pvec,tvec);
SinT=sin(T);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loop over all maps.

loop_i=0;
for i=idx_start:idx_end
  loop_i=loop_i+1;

  current_map_root=map_root_vec{i};

  br0_file_name=[current_map_root '_BR0.hdf'];
  br0_file_name_fullpath=[data_dir 'BR0/' br0_file_name];
  
  % Read in BR0 map:
  [BR0,~,~] = rdhdf_2d(br0_file_name_fullpath);
  BR0 = double(BR0);
  
  
  % Compute total signed flux in each hole using FP map for mask
  signed_flux_per_hole   = zeros(num_CHMS,1);
  area_per_hole          = zeros(num_CHMS,1);
  unsigned_flux_per_hole = zeros(num_CHMS,1);  
  
  for chm_idx=1:num_CHMS
    CH_IDX = double(CHMS==chm_idx);

    int_tmp = CH_IDX.*CHM.*BR0.*SinT.*(r_s)^2;
    int_tmp = trapz(tvec,int_tmp,1);
    signed_flux_per_hole(chm_idx) = trapz(pvec,int_tmp);   
    
    int_tmp = CH_IDX.*CHM.*abs(BR0).*SinT.*(r_s)^2;
    int_tmp = trapz(tvec,int_tmp,1);
    unsigned_flux_per_hole(chm_idx) = trapz(pvec,int_tmp);       
    
    int_tmp = CH_IDX.*CHM.*SinT.*(r_s)^2;
    int_tmp = trapz(tvec,int_tmp,1);
    area_per_hole(chm_idx) = trapz(pvec,int_tmp);   
  end  
  
  total_area_allholes          = sum(area_per_hole);
  total_open_flux_allholes     = sum(abs(signed_flux_per_hole));
  total_flux_allholes          = sum(signed_flux_per_hole);
  total_unsigned_flux_allholes = sum(unsigned_flux_per_hole);
  
  Br1au=total_open_flux_allholes/(4*pi*r1au^2);
  %Convert to nT:
  Br1auNT=Br1au*gauss2nt;
  
  disp([sprintf('%30s ',current_map_root),...
        sprintf('%12.2f ',(1e-21).*total_area_allholes),...
        sprintf('%12.2f ',(1e-22).*total_unsigned_flux_allholes),...
        sprintf('%12.2f ',(1e-22).*total_open_flux_allholes),...
        sprintf('%12.1f ',100*abs(total_flux_allholes/total_unsigned_flux_allholes)),...
        sprintf('%12.2f ',Br1auNT)...
        ]); 
end

diary off;

