function BCs  = boundary_conditions(x,g,v,u,model)


switch model.ped
    case 'pedestalmodel'
        BCs = pedestalmodel(x,g,v,u,model);
    otherwise
        [BCs.psiValue, BCs.psiRho, BCs.dpsiValue_du] = eval_BC(v,u,model.psi.BC);
        [BCs.teValue,  BCs.teRho,  BCs.dteValue_du ] = eval_BC(v,u,model.te.BC);
        [BCs.tiValue,  BCs.tiRho,  BCs.dtiValue_du ] = eval_BC(v,u,model.ti.BC);
        [BCs.neValue,  BCs.neRho,  BCs.dneValue_du ] = eval_BC(v,u,model.ne.BC);
        [BCs.niValue,  BCs.niRho,  BCs.dniValue_du ] = eval_BC(v,u,model.ni.BC);
        
        BCs.dpsiValue_dx = zeros(1,numel(x));
        BCs.dteValue_dx = zeros(1,numel(x));
        BCs.dtiValue_dx = zeros(1,numel(x));
        BCs.dneValue_dx = zeros(1,numel(x));
        BCs.dniValue_dx = zeros(1,numel(x));
end
return


function [BCs] = pedestalmodel(x,g,v,u,model)
% first test of pedestal model

[BCs.psiValue, BCs.psiRho, BCs.dpsiValue_du] = eval_BC(v,u,model.psi.BC); % psi BC not affected by pedestal mdoel

[betap,dbetap_dx] = eval_betap(x,g,v,model);
[li3,dli3_dx]  = eval_li3(x,g,v,model);
shaf = (betap+li3(end)/2);
shaf0 = 0.9375;
dshaf_dx = dbetap_dx + dli3_dx(end,:)/2;

coeffTe = 1*1600;
BCs.teValue = 1.025e3 + coeffTe*(shaf-shaf0);
BCs.teRho = 0.85;
BCs.dteValue_dx = coeffTe * dshaf_dx;
BCs.dteValue_du = zeros(1,model.dims.nu);

coeffne = 0*1e20;
BCs.neValue = 8.5e19 + coeffne*(shaf-shaf0);
BCs.neRho = 0.85;
BCs.dneValue_dx = coeffne * dshaf_dx;
BCs.dneValue_du = zeros(1,model.dims.nu);

BCs.tiValue = BCs.teValue;
BCs.tiRho   = BCs.teRho;
BCs.dtiValue_dx = BCs.dteValue_dx;
BCs.dtiValue_du = BCs.dteValue_du;

BCs.niValue  = BCs.neValue; 
BCs.niRho    = BCs.neRho;  
BCs.dniValue_dx = BCs.dneValue_dx;
BCs.dniValue_du = BCs.dneValue_du;

return