planck = 6.6256e-34;   %Plancks constant J sec
c = 2.998e8;           %speed of light m/sec
hc = 1.0/(planck*c);
oavo = 1.0/6.023e23;   %1/Avogadros number
hcoavo = hc*oavo;
RT_wbRefWLs=[400:25:700];

for il=1:13,
 if (il<10),
  fid=fopen(['ecco_oasim2_edp0',num2str(il),'_below.bin'],'r','b');
  fid2=fopen(['ecco_oasim2_esp0',num2str(il),'_below.bin'],'r','b');
 else
  fid=fopen(['ecco_oasim2_edp',num2str(il),'_below.bin'],'r','b');
  fid2=fopen(['ecco_oasim2_esp',num2str(il),'_below.bin'],'r','b');
 end
 field=fread(fid,'float32'); fclose(fid);
 field=reshape(field,360,160,12);
 field2=fread(fid2,'float32'); fclose(fid2);
 field2=reshape(field,360,160,12);
%
 RT_field_W(:,:,:,il)=field+field2;
end

for it=1:12,
  tmpW=squeeze(RT_field_W(:,:,it,:));
  for il=1:13,
   rlamm=RT_wbRefWLs(il)*1e-9;  %lambda in m
   RT_field_Eins(:,:,il) = squeeze(tmpW(:,:,il))*1e6*rlamm*hcoavo;      %Watts to uEin/s conversion
  end
  PAR_field_Eins(:,:,it)=squeeze(sum(RT_field_Eins,3));
end

fid=fopen(['ecco_oasim2_total_muEin_below.bin'],'w','b');
fwrite(fid,PAR_field_Eins,'float32');
