%INITIALIZATION class definition
%
%   Usage:
%      initialization=initialization();

classdef initialization
	properties (SetAccess=public) 
		vx            = NaN;
		vy            = NaN;
		vz            = NaN;
		vel           = NaN;
		pressure      = NaN;
		temperature   = NaN;
		waterfraction = NaN;
		sediment_head = NaN;
		epl_head      = NaN;
		epl_thickness = NaN;
		watercolumn   = NaN;
	end
	methods
        function createxml(obj,fid) % {{{
            fprintf(fid, '\n\n');
            fprintf(fid, '%s\n', '<!-- initialization -->');
            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Initial field values">','<section name="initialization" />');                    
            
			 fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="vx" type="',class(obj.vx),'" default="',obj.vx,'">','     <section name="initialization" />','     <help> x component of velocity [m/yr] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="vy" type="',class(obj.vy),'" default="',obj.vy,'">','     <section name="initialization" />','     <help> y component of velocity [m/yr] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="vz" type="',class(obj.vz),'" default="',obj.vz,'">','     <section name="initialization" />','     <help> z component of velocity [m/yr] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="vel" type="',class(obj.vel),'" default="',obj.vel,'">','     <section name="initialization" />','     <help> velocity norm [m/yr] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="pressure" type="',class(obj.pressure),'" default="',obj.pressure,'">','     <section name="initialization" />','     <help> pressure field [Pa] </help>','</parameter>');
        	 fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="temperature" type="',class(obj.temperature),'" default="',obj.temperature,'">','     <section name="initialization" />','     <help> fraction of water in the ice </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="waterfraction" type="',class(obj.waterfraction),'" default="',obj.waterfraction,'">','     <section name="initialization" />','     <help> ice thickness [m] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="sediment_head" type="',class(obj.sediment_head),'" default="',obj.sediment_head,'">','     <section name="initialization" />','     <help> sediment water head of subglacial system [m] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="epl_head" type="',class(obj.epl_head),'" default="',obj.epl_head,'">','     <section name="initialization" />','     <help> epl water head of subglacial system [m] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n','<parameter key ="watercolumn" type="',class(obj.watercolumn),'" default="',obj.watercolumn,'">','     <section name="initialization" />','     <help> thickness of subglacial water [m] </help>','</parameter>');
            fprintf(fid,'%s\n%s\n','</frame>');
        end % }}}
		function obj = initialization(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				otherwise
					error('constructor not supported');
			end
		end % }}}
		function obj = setdefaultparameters(obj) % {{{

		end % }}}
		function md = checkconsistency(obj,md,solution,analyses) % {{{
			if ismember(StressbalanceAnalysisEnum(),analyses)
				if ~(isnan(md.initialization.vx) | isnan(md.initialization.vy)),
					md = checkfield(md,'fieldname','initialization.vx','NaN',1,'size',[md.mesh.numberofvertices 1]);
					md = checkfield(md,'fieldname','initialization.vy','NaN',1,'size',[md.mesh.numberofvertices 1]);
				end
			end
			if ismember(MasstransportAnalysisEnum(),analyses),
				md = checkfield(md,'fieldname','initialization.vx','NaN',1,'size',[md.mesh.numberofvertices 1]);
				md = checkfield(md,'fieldname','initialization.vy','NaN',1,'size',[md.mesh.numberofvertices 1]);
			end
			if ismember(BalancethicknessAnalysisEnum(),analyses) & solution==BalancethicknessSolutionEnum,
				md = checkfield(md,'fieldname','initialization.vx','NaN',1,'size',[md.mesh.numberofvertices 1]);
				md = checkfield(md,'fieldname','initialization.vy','NaN',1,'size',[md.mesh.numberofvertices 1]);
				%Triangle with zero velocity
				if any(sum(abs(md.initialization.vx(md.mesh.elements)),2)==0 & sum(abs(md.initialization.vy(md.mesh.elements)),2)==0)
					md = checkmessage(md,'at least one triangle has all its vertices with a zero velocity');
				end
			end
			if ismember(ThermalAnalysisEnum(),analyses),
				md = checkfield(md,'fieldname','initialization.vx','NaN',1,'size',[md.mesh.numberofvertices 1]);
				md = checkfield(md,'fieldname','initialization.vy','NaN',1,'size',[md.mesh.numberofvertices 1]);
				if dimension(md.mesh)==3
					md = checkfield(md,'fieldname','initialization.vz','NaN',1,'size',[md.mesh.numberofvertices 1]);
				end
				md = checkfield(md,'fieldname','initialization.pressure','NaN',1,'size',[md.mesh.numberofvertices 1]);
			end
			if (ismember(EnthalpyAnalysisEnum(),analyses) & md.thermal.isenthalpy)
				md = checkfield(md,'fieldname','initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices 1]);
				md = checkfield(md,'fieldname','initialization.watercolumn'  ,'>=',0,'size',[md.mesh.numberofvertices 1]);
			end
			if ismember(HydrologyShreveAnalysisEnum(),analyses),
				if isa(md.hydrology,'hydrologyshreve'),
					md = checkfield(md,'fieldname','initialization.watercolumn','NaN',1,'size',[md.mesh.numberofvertices 1]);
				end
			end
			if ismember(HydrologyDCInefficientAnalysisEnum(),analyses),
				if isa(md.hydrology,'hydrologydc'),
					md = checkfield(md,'fieldname','initialization.sediment_head','NaN',1,'size',[md.mesh.numberofvertices 1]);
				end
			end
			if ismember(HydrologyDCEfficientAnalysisEnum(),analyses),
				if isa(md.hydrology,'hydrologydc'),
					if md.hydrology.isefficientlayer==1,
						md = checkfield(md,'fieldname', ...
														'initialization.epl_head','NaN',1,'size',[md.mesh.numberofvertices 1]);
						md = checkfield(md,'fieldname','initialization.epl_thickness','NaN',1,'size',[md.mesh.numberofvertices 1]);
					end
				end
			end
		end % }}}
		function disp(obj) % {{{
			disp(sprintf('   initial field values:'));

			fielddisplay(obj,'vx','x component of velocity [m/yr]');
			fielddisplay(obj,'vy','y component of velocity [m/yr]');
			fielddisplay(obj,'vz','z component of velocity [m/yr]');
			fielddisplay(obj,'vel','velocity norm [m/yr]');
			fielddisplay(obj,'pressure','pressure field [Pa]');
			fielddisplay(obj,'temperature','temperature [K]');
			fielddisplay(obj,'waterfraction','fraction of water in the ice');
			fielddisplay(obj,'sediment_head','sediment water head of subglacial system [m]');
			fielddisplay(obj,'epl_head','epl water head of subglacial system [m]');
			fielddisplay(obj,'epl_thickness','epl layer thickness [m]');
			fielddisplay(obj,'watercolumn','thickness of subglacial water [m]');

		end % }}}
		function marshall(obj,md,fid) % {{{

			yts=365.0*24.0*3600.0;

			WriteData(fid,'data',obj.vx,'format','DoubleMat','mattype',1,'enum',VxEnum,'scale',1./yts);
			WriteData(fid,'data',obj.vy,'format','DoubleMat','mattype',1,'enum',VyEnum,'scale',1./yts);
			WriteData(fid,'data',obj.vz,'format','DoubleMat','mattype',1,'enum',VzEnum,'scale',1./yts);
			WriteData(fid,'data',obj.pressure,'format','DoubleMat','mattype',1,'enum',PressureEnum);
			WriteData(fid,'data',obj.temperature,'format','DoubleMat','mattype',1,'enum',TemperatureEnum);
			WriteData(fid,'data',obj.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum);
			WriteData(fid,'data',obj.sediment_head,'format','DoubleMat','mattype',1,'enum',SedimentHeadEnum);
			WriteData(fid,'data',obj.epl_head,'format','DoubleMat','mattype',1,'enum',EplHeadEnum);
			WriteData(fid,'data',obj.epl_thickness,'format','DoubleMat','mattype',1,'enum',HydrologydcEplThicknessEnum);
			WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum);

			if md.thermal.isenthalpy,
				tpmp = md.materials.meltingpoint - md.materials.beta*md.initialization.pressure;
				pos  = find(md.initialization.temperature>tpmp);
				enthalpy      = md.materials.heatcapacity*(md.initialization.temperature-md.constants.referencetemperature);
				enthalpy(pos) = md.materials.heatcapacity*tpmp(pos) - md.constants.referencetemperature + md.materials.latentheat*md.initialization.waterfraction(pos);
				WriteData(fid,'data',enthalpy,'format','DoubleMat','mattype',1,'enum',EnthalpyEnum());
			end
		end % }}}
	end
end
