%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This program is used to run the GDL model 
% Jason Siegel, Denise McKay
% Summer 2007
%
% Loads experimental data and runs the model for a
% given [alpha_dw, t_wl] pair, saves the model inputs and outputs to a .mat
% file to be used by Voltage_parameterization.m to tune the voltage model
%
% This program calls the following files:
% fcs_params.m                  ~ parameter values
% FileName.mat                  ~ calibration experimental data
% mscaledata.m                  ~ scales experimental data into SI units
% mscaledata_a.m                ~ scales experimental data into non-SI units
% satpressure.m                 ~ given temperature, finds water vapor
%                                 sat pressure
% GDLModel.sml                  ~ GDL model built in Simulink
% GDL_Mass_balance.m            ~ checks mass balance in each volume
%                                 (optional)
% plot_anode.m                  ~ plots anode variables
%                                 (optional)
% plot_cathode.m                ~ plots cathode variables
%                                 (optional)
% plot_membrane.m               ~ plots membrane variables
%                                 (optional)
% plot_model_inputs.m           ~ plots experimental data
%                                 (optional)
% plot_outputs.m                ~ plots model output variables
%                                 (optional)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
close all;
clc
addpath MatFiles

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nomenclature
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% t         ~ simulation time (sec)
% x         ~ simulated states
% y         ~ simulated outputs, used for voltage param ID
%               AveV,Tst,c_o2,c_h2,p_ca         %1,2,3,4,5
%               lambda,I_st,m_l_an(4),s_ca(1)   %6,7,8,9               
%
% time      ~ experimental time (sec)
% ur        ~ formed from mscaledata, removed time vector
%               PCaIn,TCaIn,RHCaIn,PAnIn,TAnIn,RHAnIn       %1,2,3,4,5,6
%               PCaOut,TCaOut,RHCaOut,PAnOut,TAnOut,RHAnOut %7,8,9,10,11,12
%               WaCaIn,IFC,TFC,PurgeState                   %13,14,15,16 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load Experimental Data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Can specify either a .xls (load_xls=1) or .mat (load_xls=0) file format

load_xls=0;     % use file format other than .xls 
%load_xls=1;    % use .xls file format 

if(load_xls)
    FileName='12.21.2004_11_09.xls';
    skip_start=9;skip_stop=0; % number of data rows removed from begining and end of file
    data=open_file(FileName,skip_start,skip_stop); % open excel spreadsheet, clip data
else
    FileName='Dec18FulldayData';
    load(FileName,'data');
end

% scale data into appropriate units
ExperData_set=mscaledata(data); % scales experimental data to desired units
calcdata % run script to calculate data based on measurements
[nrows,ncols]=size(ExperData_set);
ur = ExperData_set(:,2:ncols);
time = ExperData_set(:,1); % experimental time (sec)
clear ExperData_set

% load system parameters
fcs_params

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specify Initial Conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% measurements used for initialization
T=ur(1,15);         % cell temperature (K)
PCaIn_Run=ur(1,1);  % cathode inlet pressure (Pa)
PCaOut_Run=ur(1,7); % cathode outlet pressure (Pa)
PAnIn_Run=ur(1,4);  % anode inlet pressure

% initial channels concentrations, specify using ideal gas law
cv_IC=0.99*satpressure(T)/(R_G_U*T);
cH2IC=(PAnIn_Run-satpressure(T))/(R_G_U*T);
cO2IC=((PCaIn_Run+PCaOut_Run)/2-satpressure(T))/(R_G_U*T)*OMF_ca_in;
cN2IC=((PCaIn_Run+PCaOut_Run)/2-satpressure(T))/(R_G_U*T)*(1-OMF_ca_in);

% calculate initial channels masses given initial concentrations
mH2IC=cH2IC*FC_ANODE_VOLUME*HYDROGENMOLARMASS;
mH2IC_OM=cH2IC*FC_ANODE_OM_VOLUME*HYDROGENMOLARMASS;
mWIC_AN=cv_IC*FC_ANODE_VOLUME*VAPORMOLARMASS;
mWIC_AN_OM=cv_IC*FC_ANODE_OM_VOLUME*VAPORMOLARMASS;
mWIC_CA=cv_IC*FC_CATHODE_VOLUME*VAPORMOLARMASS;
mO2IC=cO2IC*FC_CATHODE_VOLUME*OXYGENMOLARMASS;
mN2IC=cN2IC*FC_CATHODE_VOLUME*NITROGENMOLARMASS;

% initial s in each GDL section, set at immobile saturation
sIC=imobile_sat;
sIC1=imobile_sat;
sIC2=imobile_sat;
sIC3=imobile_sat;

% initial anode GDL masses given s in each GDL section
mH2IC_3=cH2IC*(1-sIC3)*V_P;
mv_IC_3=cv_IC*(1-sIC3)*V_P;
mH2IC_2=cH2IC*(1-sIC2)*V_P;
mv_IC_2=cv_IC*(1-sIC2)*V_P;
mH2IC_1=cH2IC*(1-sIC1)*V_P;
mv_IC_1=cv_IC*(1-sIC1)*V_P;

% initial cathode GDL masses given s in each GDL section
% Note, there is no N2 as c_N2 is constant and not a state
mO2IC=cO2IC*(1-sIC)*V_P;
mv_IC=cv_IC*(1-sIC)*V_P;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run the GDL model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model_to_run='GDLModel';

% set solver criteria
optimset = simset('solver','ode23s','RelTol','auto','AbsTol','auto');%,'MinStep',1e-13);
deltat=-1;      % simulation sample time step

% to simulate only a portion of the calibration data set
if(strcmp(FileName,'Dec18FulldayData'))
    indexx=find(time<15200); % use only part of this data set for simulation
else
    indexx=1:length(time); % use entire data set for simulation
end
Ldata=max(indexx); % index of last data point to use in simulation

% run simulation for specified time period
[t,x,y] = sim(model_to_run,[time(indexx).'],optimset,[time,ur]);

% % perform mass balance on each volume
%GDL_Mass_balance(tout,GDL_out)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Specify Model Outputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Transport Diagram  ("g" used to denote gas)
%                 _________________________________________
%                |   GDL 1     |   GDL 2     |   GDL 3     |  
%                |  s_e(1)     |   s_e(2)    |   s_e(3)    
%   Membrane     |  c_g_ca(1)  |   c_g_ca(2) |   c_g_ca(3) |      Channel
%                |  cv_e(1)    |   cv_e(2)   |   cv_e(3)   |
%                |             |             |             |
%                |          W_g_e(1)-->   W_g_e(2)-->   W_g_e(3)--> 
%                |          Wv_e(1)-->     Wv_e(2)-->     Wv_e(3)-->
%                |             |             |             |
%                |         W_l_e(1)-->     W_l_e(2)-->    W_l_e(3)-->
%                |_____________|_____________|_____________|
%                -------> y-direction
%
% Order of variables in CATHODE_k           Order of variables in ANODE_k 
% 1. c_o2(k)                                1. c_h2(k)
% 2. c_n2(k)                                2. c_v(k)
% 3. c_v(k)                                 3. W_h2 (k)
% 4. W_o2 (k)                               4. W_v(k)
% 5. W_n2 (k)                               5. s_an(k)
% 6. W_v(k)                                 6. W_l_an(k)
% 7. s_ca(k)                                7. T_an(k)
% 8. W_l_ca(k)
% 9. T_ca(k)

% Order of variables in MEMBR_out
% 1. lambda_mb
% 2. lambda_an
% 3. lambda_ca
% 4. WH2_react
% 5. WO2_react
% 6. Wv_gen
% 7. Wv_mb
% 8. Nv_drag
% 9. Nv_diff

% Order of variables in ANODE_Pressure     Variables in CATHODE_Pressure
% 1. p_an,ch                               1. p_ca,ch
% 2. p_h2,ch                               2. p_v,ca,ch
% 3. p_v,an,ch                             3. p_o2,ca,ch
% 4. p_an,om                               4. p_n2,ca,ch

% Order of variables in Current_Density
% 1. i_app
% 2. i

% Order of variables in Excess_Ratio
% 1. W_H2_IN
% 2. W_H2_rct (total stack)
% 3. W_O2_IN
% 4. W_O2_rct (total stack)

% oxygen partial pressures  (bar)
p_o2=[];
p_o2(:,1)=CATHODE_1(:,1).*R_G_U.*CATHODE_1(:,9)/100000;  
p_o2(:,2)=CATHODE_2(:,1).*R_G_U.*CATHODE_2(:,9)/100000;  
p_o2(:,3)=CATHODE_3(:,1).*R_G_U.*CATHODE_3(:,9)/100000;  
p_o2(:,4)=CATHODE_Pressure(:,3)/100000;

% oxygen partial pressures  (bar)
p_n2=[];
p_n2(:,1)=CATHODE_1(:,2).*R_G_U.*CATHODE_1(:,9)/100000;  
p_n2(:,2)=CATHODE_2(:,2).*R_G_U.*CATHODE_2(:,9)/100000;  
p_n2(:,3)=CATHODE_3(:,2).*R_G_U.*CATHODE_3(:,9)/100000;  
p_n2(:,4)=CATHODE_Pressure(:,4)/100000;

% cathode water vapor partial pressures  (bar)
p_v_ca=[];
p_v_ca(:,1)=CATHODE_1(:,3).*R_G_U.*CATHODE_1(:,9)/100000;
p_v_ca(:,2)=CATHODE_2(:,3).*R_G_U.*CATHODE_2(:,9)/100000;
p_v_ca(:,3)=CATHODE_3(:,3).*R_G_U.*CATHODE_3(:,9)/100000;
p_v_ca(:,4)=CATHODE_Pressure(:,2)/100000;

% hydrogen partial pressures  (bar)
p_h2=[];
p_h2(:,1)=ANODE_1(:,1).*R_G_U.*ANODE_1(:,7)/100000;  
p_h2(:,2)=ANODE_2(:,1).*R_G_U.*ANODE_2(:,7)/100000; 
p_h2(:,3)=ANODE_3(:,1).*R_G_U.*ANODE_3(:,7)/100000;
p_h2(:,4)=ANODE_Pressure(:,2)/100000; 

% anode water vapor partial pressures (bar)
p_v_an=[];
p_v_an(:,1)=ANODE_1(:,2).*R_G_U.*ANODE_1(:,7)/100000;
p_v_an(:,2)=ANODE_2(:,2).*R_G_U.*ANODE_2(:,7)/100000;
p_v_an(:,3)=ANODE_3(:,2).*R_G_U.*ANODE_3(:,7)/100000;
p_v_an(:,4)=ANODE_Pressure(:,3)/100000;

% anode relative humidity, not constrained to 1
RH_an=[];
RH_an(:,1)=p_v_an(:,1)./(satpressure(ANODE_1(:,7))/100000);
RH_an(:,2)=p_v_an(:,2)./(satpressure(ANODE_2(:,7))/100000);
RH_an(:,3)=p_v_an(:,3)./(satpressure(ANODE_3(:,7))/100000);

% cathode relative humidity, not constrained to 1
RH_ca=[];
RH_ca(:,1)=p_v_ca(:,1)./(satpressure(CATHODE_1(:,9))/100000);
RH_ca(:,2)=p_v_ca(:,2)./(satpressure(CATHODE_2(:,9))/100000);
RH_ca(:,3)=p_v_ca(:,3)./(satpressure(CATHODE_3(:,9))/100000);

% anode mass of liquid water in anode GDL sections (g)
m_l_an=[];
m_l_an(:,1)=ANODE_1(:,5)*WATERDENSITY*V_P*1000;
m_l_an(:,2)=ANODE_2(:,5)*WATERDENSITY*V_P*1000;
m_l_an(:,3)=ANODE_3(:,5)*WATERDENSITY*V_P*1000;
m_l_an(:,4)=y(:,8)*1000;

% anode mass of liquid water in cathode GDL sections (g)
m_l_ca=[];
m_l_ca(:,1)=CATHODE_1(:,7)*WATERDENSITY*V_P*1000;
m_l_ca(:,2)=CATHODE_2(:,7)*WATERDENSITY*V_P*1000;
m_l_ca(:,3)=CATHODE_3(:,7)*WATERDENSITY*V_P*1000;

% voltage
V_model =y(:,1); % predicted avg cell voltage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mscaledata_a;       % scales experimental data into non-SI units
%plot_anode          % plot anode variables
%plot_cathode        % plot cathode variables
%plot_membrane       % plot membrane variables
plot_outputs        % plot model output variables

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Save Data to File for Offline Parameter ID
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% save output matrix y to a .mat file
%save Dec18_adw12_0_tlw1200 % enter name of .mat file