close all;
clear all;

r_s=6.957e10;

frmtstr='%11.3f ';

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

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' ...
};

output_file_name='ofp_table1_fluxes_polefields_from_BR0';

idx_start=1;
idx_end=length(map_root_vec);

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

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 ','MAP NAME'),...
             sprintf('%12s ','UnsignedFlux')...
             sprintf('%12s ','BrA65(G)'),...
             sprintf('%12s ','BrB65(G)'),...
        ]);  
  end

  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,tvec,pvec] = rdhdf_2d(br0_file_name_fullpath);
  BR0 = double(BR0);
  tvec = double(tvec);
  pvec = double(pvec);

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

  % Compute integrand:
  br_int=BR0.*SinT*(r_s)^2;
  
  % Signed flux:
  flux_signed=trapz(tvec,br_int,1);
  flux_signed=trapz(pvec,flux_signed);
  
  % Unsigned flux:
  flux_unsigned=trapz(tvec,abs(br_int),1);
  flux_unsigned=trapz(pvec,flux_unsigned);

  % Positive flux:
  br_pos=br_int;
  br_pos(br_pos<0)=0;
  flux_positive=trapz(tvec,br_pos,1);
  flux_positive=trapz(pvec,flux_positive);

  % Negative flux:
  br_neg=br_int;
  br_neg(br_neg>0)=0;
  flux_negative=trapz(tvec,br_neg,1);
  flux_negative=trapz(pvec,flux_negative);

  %Average mag field above 65 degrees:
  area_int=ones(size(br_int)).*SinT*(r_s)^2;
  area_int(T<(90+65)*pi/180)=0;
  area_above_65=trapz(tvec,area_int,1);
  area_above_65=trapz(pvec,area_above_65);
  br_pos=br_int;
  br_pos(T<(90+65)*pi/180)=0;
  br_above65=trapz(tvec,br_pos,1);
  br_above65=trapz(pvec,br_above65);  
  br_above65=br_above65./area_above_65;

  %Average mag field below 65 degrees:
  area_int=ones(size(br_int)).*SinT*(r_s)^2;
  area_int(T>25*pi/180)=0;
  area_below_65=trapz(tvec,area_int,1);
  area_below_65=trapz(pvec,area_below_65);
  br_pos=br_int;
  br_pos(T>25*pi/180)=0;
  br_below65=trapz(tvec,br_pos,1);
  br_below65=trapz(pvec,br_below65);  
  br_below65=br_below65./area_below_65;    
    
  disp([sprintf('%50s ',br0_file_name),...
        sprintf('%12.1f ',1e-22*flux_unsigned),...
        sprintf('%12.1f ',br_above65),...
        sprintf('%12.1f ',br_below65)...
        ]);      

end
    
%Close output text file:
diary off;

