close all;
clear all;

%Mx=Gauss-cm^2
gauss2nt=1e5;
r_s=6.957e10;
r1au=1.496e13;

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

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_fluxes_IMF_from_BR1';

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 ','Br1AU(nT)')... 
        ]);  
  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 R values
  for ir1=1:length(r1_vec)
      
    r1=r1_vec(ir1);

    if (r1<30)
      br1_file_name=[current_map_root '_PFSS' num2str(r1*10,'%2d') '_BR1.hdf'];
    else
      br1_file_name=[current_map_root '_MHD_BR1.hdf'];   
    end
    br1_file_name_fullpath=[data_dir 'BR1/' br1_file_name];
  
    % Read in BR1 map:
    [BR1,br_t,br_p] = rdhdf_2d(br1_file_name_fullpath);
    BR1 = double(BR1);
    br_t = double(br_t);
    br_p = double(br_p);

    % Get 2D scale information and sint:
    [P,T]=meshgrid(br_p,br_t);
    SinT=sin(T);

    % Compute integrand:
    br_int=BR1.*SinT*(r_s*r1)^2;
  
    % Signed flux:
    flux_signed=trapz(br_t,br_int,1);
    flux_signed=trapz(br_p,flux_signed);
  
    % Unsigned flux:
    flux_unsigned=trapz(br_t,abs(br_int),1);
    flux_unsigned=trapz(br_p,flux_unsigned);

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

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

    %Br value at 1AU based on unsigned flux:
    %Get Br at 1AU (gauss)
    Br1au=flux_unsigned./(4*pi*r1au^2);
    %Convert to nT:
    Br1au=Br1au*gauss2nt;
  
    disp([sprintf('%50s ',br1_file_name),...
          sprintf('%12.2f ',Br1au)...
          ]);
  end
end
    
%Close output text file:
diary off;

