function [f] = fsagsesolv(g2,g2prime,g3,g3prime, iota, iotaprime, pprime, RBphivac, Lam,Lamp, solveopt)
% Solve flux surface averaged grad-shafranov equation for F
%
% INPUTS
% g2,g2prime: values and psi derivative of <gradV^2/R^2>
% g3,g3prime: values and psi derivative of <1/R^2>
% q, qprime : values and psi derivative of q
% pprime    : value of dp/dpsi
% RBphivac  : vacuum R*Bphi;
% Lam, Lamp : Finite element matrices (and first spatial derivatives) such
%               that F.^2 = Lam*f, 2FF' = Lamp*f;
% solveopt  : 'H2' to solve for H^2 (H=dpsi/dV), 'F2' to solve for F^2 (F=RBphi);
% 
% OUTPUTS: f: Finite element coefficients for F^2, such that
%                   F = sqrt(Lam*f); 
%                   FF' = 1/2*Lamp*f;
%                   or similar on other grids

switch solveopt
    case 'F2'
        a = 1/2*((g3.*iota).^2.*g2 + 16*pi^4*g3);
        b = g2prime .* (g3.*iota).^2 + g2.*(g3.*iota) .* (g3.*iotaprime + g3prime.*iota);
        c = 16*pi^4*4e-7*pi * pprime; 
    case 'H2'
        a = 1/2*(g2 + 16*pi^4 * 1./(g3.*iota.^2) );
        b = 16*pi^4 * 1./iota .* (-1./((g3.*iota).^2) .* (g3prime.*iota + g3.*iotaprime))  + g2prime;
        c = 4*pi.^2 * 4e-7*pi .* pprime;      
end

B = -b ./ (a);
C = -c ./ (a);

% solve first order linear ODE
% where y = F.^2 or H.^2;
% y'(x) = C(x) + B(x)y(x)

%% Finite element solution - formulate as algebric matrix-vector equation

% solve ODE with constraint on edge
A = (-bsxfun(@times,B,Lam) + Lamp);
b = C;

% constraint Lam(end,:)*f = bc
% null of constraint
f(size(Lam,2),:) = RBphivac.^2; % last one fixed by boundary condition
f(1:end-1,:) = A(:,1:end-1)\(b-A(:,end)*f(end)); % solve the rest

check = 0;
if check
    % check solution
    F2 = Lam*f;
    rhsF2 = Lamp*f;
    lhsF2 = C + B.*F2;
    figure(55); clf;

    subplot(211);
    plot(F2,'k-'); hold on; grid on; 
    plot([0 1],RBphivac.^2*[1 1],'ok');

    subplot(212);
    plot(rhsF2,'ko'); hold on;
    plot(lhsF2,'kx'); grid on;
    
    drawnow; pause(0.05);
end

%%
return