disp('this tutorial is still under construction');
return % return so that testing scripts don't error on this

%%
close all hidden;
run(fullfile(pwd,'..','RAPTOR_path.m')); % add RAPTOR path
%%
[config] = RAPTOR_config; % load default params
%%
[model,params,init,g,v,U] = build_RAPTOR_model(config);

%%
x0 = RAPTOR_initial_conditions(model,init,g(:,1),v(:,1)); 

%% 
% Now we are ready to run RAPTOR - this will display some plasma quantities as
% they evolve.
simres = RAPTOR_predictive(x0,g,v,U,model,params);

% To compute some physical outputs, we run the RAPTOR_out function.
out = RAPTOR_out(simres,model,params);


%% Find stationary condition
x0 = simres.X(:,end);
u0 = simres.U(:,end);
xdot0 = zeros(size(x0));
uplgauss = ones(model.rgrid.nrhogauss,1)*out.upl(end,end);
xdot0(model.psi.xind,:) = model.psi.Lamgauss\uplgauss;

x0mask = false(size(x0)); x0mask(end) = true; % fix last element of x
u0mask = true(size(u0)); u0mask(1,:) = false; % free Ip
xdot0mask = true(size(x0));
fixmask= {x0mask,u0mask,xdot0mask};

[x,u,xdot,df_dxdot,df_dx,df_du,f,stap,geop,trap] = RAPTOR_solve_state(x0,u0,xdot0,g,v,model,params,fixmask);
%%
out = RAPTOR_out(x,g,v,xdot,u,1,stap,geop,trap,model,params);
plot(out.rho,out.q)

