close all;
clear all;

r_s=6.957e10;

frmtstr='%11.3f ';

data_dir = '../data/';
results_dir = '../results/';

chm_S_obs_file_name      = [data_dir 'CHM/' 'CR2098_S_20100708_OBS_CHM.hdf'];
no_data_region_file_name = [data_dir 'CHM/' 'CR2098_S_20100708_OBS_NODATAMASK.hdf'];

map_root_vec={ ...
'CR2098_S_20100707_ADAPT_FAR' ...
'CR2098_S_20100707_ADAPT_FLIPFAR' ...
'CR2098_S_20100707_ADAPT_NOFAR' ...
'CR2098_S_20100708_GONG'...
'CR2098_D_GONG' ...
'CR2098_S_20100708_HMI_LOS' ...
'CR2098_D_HMI_LOS' ...
'CR2098_D_HMI_VECTOR' ...
'CR2098_S_20100708_LMSALESFAM' ...
'CR2098_S_20100708_MDI' ...
'CR2098_D_MDI' ...
'CR2098_D_VSM' ...
'CR2098_D_VSM_PFIT' ...
'CR2098_S_20100708_VSM' ...
};

mhd_runs_idx=[1 2 3 9 13];   %Indices of above maps where there is MHD data.
pfss_scaling_idx=12;         %Index of run with extra pfss radii.

output_file_name='ofp_table1_obsSchm_vs_ofm';

idx_start=1;
idx_end=length(map_root_vec);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

% Load observed CHM NODATA mask:
[ND,~,~] = rdhdf_2d(no_data_region_file_name);
ND = double(ND);  
NDMASK=abs(ND-1);
  
% Get 2D scale information and sint:
tvec = double(tvec);
pvec = double(pvec);
[~,T]=meshgrid(pvec,tvec);
SinT=sin(T);

% Get area of observed holes (automatically excludes no-data region).
CHM_AREA_INT=CHM.*SinT.*(r_s)^2;
CHM_AREA=trapz(tvec,CHM_AREA_INT,1);
CHM_AREA=trapz(pvec,CHM_AREA);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loop over all maps and their runs.

loop_i=0;
for i=idx_start:idx_end
  loop_i=loop_i+1;
  
  % Output data to text file:
  if(loop_i==1)
    [~,~]=system(['rm ' results_dir output_file_name '.txt']);
    diary([results_dir output_file_name '.txt']);
    disp([sprintf('%50s ','OBSERVED AREA:'),sprintf(frmtstr,(1e-21).*CHM_AREA)]);
    disp([sprintf('%50s ','MAP NAME'),...
      sprintf('%14s ','OF_AreaMasked'),...
      sprintf('%14s ','Amask-Aobs'),...    
     ]);  
  end

  current_map_root=map_root_vec{i};

  %If i is in mhd list, r1vec has 30 in it, otherwise not.
  if (min(abs(i-mhd_runs_idx))==0)
    r1_vec=[2.5 2.0 30];
  else
    r1_vec=[2.5 2.0];
  end
  
  if (i==pfss_scaling_idx)
    r1_vec = [r1_vec 1.4 1.3];
  end

  %LOOP OVER RUNS
  for ir1=1:length(r1_vec)

    r1=r1_vec(ir1);

    if (r1<30)
      ofm_file_name=[current_map_root '_PFSS' num2str(r1*10,'%2d') '_OFM.hdf'];
    else
      ofm_file_name=[current_map_root '_MHD_OFM.hdf'];   
    end
    ofm_file_name_fullpath=[data_dir 'OFM/' ofm_file_name];

    % Read in OFM map.
    [OFM,~,~] = rdhdf_2d(ofm_file_name_fullpath);
    OFM = double(OFM);  
    
    % Compute total area of OFM:
    OFM_AREA_INT = OFM.*SinT.*(r_s)^2;
    OFM_AREA     = trapz(tvec,OFM_AREA_INT,1);
    OFM_AREA     = trapz(pvec,OFM_AREA);
    
    % Compute no-data masked area of OFM:
    OFM_AREA_MASKED_INT = NDMASK.*OFM.*SinT.*(r_s)^2;
    OFM_AREA_MASKED     = trapz(tvec,OFM_AREA_MASKED_INT,1);
    OFM_AREA_MASKED     = trapz(pvec,OFM_AREA_MASKED);
   
    disp([sprintf('%50s ',ofm_file_name),...
          sprintf('%14.1f ',(1e-21).*OFM_AREA_MASKED),...
          sprintf('%14.1f ',(1e-21).*(OFM_AREA_MASKED-CHM_AREA)),...      
        ]);         
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

diary off;










