function [Tout,denout] =  ETM_model(Year,DoY,UTsec,Alt,Lon,Lat,F107,F107A,Aplist)


%    Year:  Year
%    DoY:  Date number of year
%    UTsec:    seconds in day in universal time (UT)
%    Alt:     altitude in km
%    Lon:   geodetic longitude in degrees.
%    Lat:    geodetic latitude in degrees.
%    F107:      daily F10.7 flux for previous day
%    F107A:   81 day average of F10.7 flux (centered on doy).
%    Aplsit:    magnetic index information
%                Aplist = [Ap1,Ap2,Ap3,Ap4]
%                                  Ap1:       daily averaged Ap
%                                  Ap2:       3 hour AP for current time
%                                  Ap3:       3 hour AP for 3 hours before current time
%                                  Ap4:       3 hour AP for 6 hours before current time


              currentFolder = pwd;
              addpath([genpath(currentFolder),'codedata'])
              
              
              
               Am = matrix_cal(F107A,F107,Aplist,DoY,UTsec);
               P     = PC_cal(UTsec,Lon,Lat);
               
               T         =  sum(Am.*P,1);                                                                                                                       % Compute exospehric temperature utilizing 5 orders of PCs 
               Tout   =  reshape(T,size(Year));
               
               den = ones(size(T));
               
               for ik = 1:numel(Year)
                     den(ik) =  nrlMSIS(Year(ik),DoY(ik),UTsec(ik),Alt(ik),Lat(ik),Lon(ik),F107A(ik),F107(ik),Aplist(ik,2),T(ik));                      % Compute the mass density with the aid of MSIS
               end
               
               denout  = reshape(den,size(Year));
end
               



