%% plasieqf.m calls plasieq2.m to determine the electron density as a function of L for
%% different MLT for the request UT. It then saves the results in the ASCII file
%% GPID_eq_YYYY_DDD_TT.txt where YYYY is the year, DDD is the day number, and TT is
%% UT to the nearest hour.
%% 
%% Author: Phillip A. Webb
%% Last updated: January 10 2005
%%

function plasieqf(timmat,dno,year,mltv)

% timmat = UT of points of interest (scalar)
% dno = day number of points of interest (scalar)
% year = year of points of interest (scalar)
% mlt = MLT of points of interest (scalar or vector)

%% Determine the correct day number string
if dno<10
 dnos=['00',int2str(dno)];
elseif dno<100
 dnos=['0',int2str(dno)];
else
 dnos=int2str(dno);
end

%% Determine the correct hour string
if timmat<10
 tims=['0',int2str(round(timmat))];
else
 tims=int2str(round(timmat));
end

% Conduct initial call to determine the L-shell values
[eden,oxden,hyden,heden,hytot,hetot,lshv]=plasieq2(mltv(1),timmat,dno,year);

%Open up the required file
filse=['GPID_',int2str(year),'_',dnos,'_',tims,'_eq.txt'];
fid=fopen(filse,'w');
tits=['GPID Electron Equatorial Density (m^-3)  UT=',tims,'  Day number=',int2str(dno),'  Year=',int2str(year)];
fprintf(fid,[tits,'\n\n']);
fprintf(fid,'Number of MLT values\n');
fprintf(fid,'%3.0f\n\n',length(mltv));
fprintf(fid,'Number of L values\n');
fprintf(fid,'%3.0f\n\n',length(lshv));
fprintf(fid,'MLT (fractional hours)\n');
fprintf(fid,'%5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n',mltv(:)');
fprintf(fid,'\n\n');
fprintf(fid,'L\n');
fprintf(fid,'%5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n',lshv(:)');
fprintf(fid,'\n\n');
for i=1:length(mltv)
 if i>1
  [eden,oxden,hyden,heden,hytot,hetot,lshv]=plasieq2(mltv(i),timmat,dno,year);
 end 
 fprintf(fid,'MLT = %5.2f \n',mltv(i));
 fprintf(fid,'%3.2e %3.2e %3.2e %3.2e %3.2e %3.2e %3.2e %3.2e\n',real(eden));
 fprintf(fid,'\n\n');
end 
fclose(fid);