% computes the sw radiative fluxes at the surface given an input absorption and reflection percent in the atmosphere
% then outpus the difference between the simple theory and the onservations--- for optimization 
function F = myfun_solar_full(x)

%x= input [reflection percent; absorption percent]

if length(x) > 2 | length(x) < 2
    error('input vector must be length 2')
end

R=x(1); %reflection percent
A=x(2); %absorption percent

load ('local_Solar.mat', 'S_down_toa', 'S_up_toa', 'S_down_surf', 'S_up_surf','absorb'); %should have S_bar, and alb_bar, CO2, RH preloaded

alb = S_up_surf/S_down_surf; %surface albedo -higher order reflections should obey the same albedo so should be exact
S_out = [ S_down_toa*R  + alb*S_down_toa*(1-R-A)^2 + R*(alb^2)*S_down_toa*(1-R-A)*(1-R-A) + (R^2)*(alb^3)*S_down_toa*(1-R-A)*(1-R-A);   %predicted upwelling at TOA  
          S_down_toa*A + A*alb*S_down_toa*(1-R-A) + A*R*(alb^2)*S_down_toa*(1-R-A) + A*(R^2)*(alb^3)*S_down_toa*(1-R-A);    %predicted absorption
          S_down_toa*(1-R-A)+R*alb*S_down_toa*(1-R-A)+(R^2)*(alb^2)*S_down_toa*(1-R-A)+(R^3)*(alb^3)*S_down_toa*(1-R-A)] ;%surface predicted downwelling
      
S_in = [S_up_toa; absorb; S_down_surf];

F=S_out-S_in;
          