%MATICE class definition
%
%   Usage:
%      matice=matice();

classdef matice
	properties (SetAccess=public) 
		rho_ice                    = 0.;
		rho_water                  = 0.;
		rho_freshwater             = 0.;
		mu_water                   = 0.;
		heatcapacity               = 0.;
		latentheat                 = 0.;
		thermalconductivity        = 0.;
		temperateiceconductivity   = 0.;
		meltingpoint               = 0.;
		beta                       = 0.;
		mixed_layer_capacity       = 0.;
		thermal_exchange_velocity  = 0.;
		rheology_B   = NaN;
		rheology_n   = NaN;
		rheology_law = '';

		%gia: 
		lithosphere_shear_modulus  = 0.;
		lithosphere_density        = 0.;
		mantle_shear_modulus       = 0.;
		mantle_density             = 0.;

	end
	methods
        function createxml(obj,fid) % {{{
            fprintf(fid, '\n\n');
            fprintf(fid, '<!-- materials -->\n');
			 fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="rho_ice" type="',class(obj.rho_ice),'" default="',convert2str(obj.rho_ice),'">','     <section name="materials" />','     <help> ice density [kg/m^3] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="rho_water" type="',class(obj.rho_water),'" default="',convert2str(obj.rho_water),'">','     <section name="materials" />','     <help> ocean water density [kg/m^3] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="rho_freshwater" type="',class(obj.rho_freshwater),'" default="',convert2str(obj.rho_freshwater),'">','     <section name="materials" />','     <help> fresh water density [kg/m^3] </help>','</parameter>');
             
  
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="mu_water" type="',class(obj.mu_water),'" default="',convert2str(obj.mu_water),'">','     <section name="materials" />','     <help> water viscosity [N s/m^2] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="heatcapacity" type="',class(obj.heatcapacity),'" default="',convert2str(obj.heatcapacity),'">','     <section name="materials" />','     <help> heat capacity [J/kg/K] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="latentheat" type="',class(obj.latentheat),'" default="',convert2str(obj.latentheat),'">','     <section name="materials" />','     <help> latent heat of fusion [J/kg] </help>','</parameter>');
             
             
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="thermalconductivity" type="',class(obj.thermalconductivity),'" default="',convert2str(obj.thermalconductivity),'">','     <section name="materials" />','     <help> ice thermal conductivity [W/m/K] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="temperateiceconductivity" type="',class(obj.temperateiceconductivity),'" default="',convert2str(obj.temperateiceconductivity),'">','     <section name="materials" />','     <help> temperate ice thermal conductivity [W/m/K] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="meltingpoint" type="',class(obj.meltingpoint),'" default="',convert2str(obj.meltingpoint),'">','     <section name="materials" />','     <help> melting point of ice at 1atm in K </help>','</parameter>');
             
             
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="beta" type="',class(obj.beta),'" default="',convert2str(obj.beta),'">','     <section name="materials" />','     <help> rate of change of melting point with pressure [K/Pa] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="mixed_layer_capacity" type="',class(obj.mixed_layer_capacity),'" default="',convert2str(obj.mixed_layer_capacity),'">','     <section name="materials" />','     <help> mixed layer capacity [W/kg/K] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="thermal_exchange_velocity" type="',class(obj.thermal_exchange_velocity),'" default="',convert2str(obj.thermal_exchange_velocity),'">','     <section name="materials" />','     <help> thermal exchange velocity [m/s] </help>','</parameter>');
             
             
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="rheology_B" type="',class(obj.rheology_B),'" default="',convert2str(obj.rheology_B),'">','     <section name="materials" />','     <help> flow law parameter [Pa/s^(1/n)] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="rheology_n" type="',class(obj.rheology_n),'" default="',convert2str(obj.rheology_n),'">','     <section name="materials" />','     <help> Glens flow law exponent </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="rheology_law" type="',class(obj.rheology_law),'" default="',convert2str(obj.rheology_law),'">','     <section name="materials" />','     <help> law for the temperature dependance of the rheology: "None", "Paterson",  "Arrhenius" or "LliboutryDuval" </help>','</parameter>');
             
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="lithosphere_shear_modulus" type="',class(obj.lithosphere_shear_modulus),'" default="',convert2str(obj.lithosphere_shear_modulus),'">','     <section name="materials" />','     <help> Lithosphere shear modulus [Pa] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="lithosphere_density" type="',class(obj.lithosphere_density),'" default="',convert2str(obj.lithosphere_density),'">','     <section name="materials" />','     <help> Lithosphere density [g/cm^-3] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="mantle_shear_modulus" type="',class(obj.mantle_shear_modulus),'" default="',convert2str(obj.mantle_shear_modulus),'">','     <section name="materials" />','     <help> Mantle shear modulus [Pa] </help>','</parameter>');
             fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="mantle_density" type="',class(obj.mantle_density),'" default="',convert2str(obj.mantle_density),'">','     <section name="materials" />','     <help> Mantle density [g/cm^-3] </help>','</parameter>');
        
         
        end % }}}
		function obj = matice(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				case 1
					inputstruct=varargin{1};
					list1 = properties('matice');
					list2 = fieldnames(inputstruct);
					for i=1:length(list1)
						fieldname = list1{i};
						if ismember(fieldname,list2),
							obj.(fieldname) = inputstruct.(fieldname);
						end
					end
				otherwise
					error('constructor not supported');
			end
		end % }}}
		function obj = setdefaultparameters(obj) % {{{

			%ice density (kg/m^3)
			obj.rho_ice=917.;

			%ocean water density (kg/m^3)
			obj.rho_water=1023.;

			%fresh water density (kg/m^3)
			obj.rho_freshwater=1000.;

			%water viscosity (N.s/m^2)
			obj.mu_water=0.001787;  

			%ice heat capacity cp (J/kg/K)
			obj.heatcapacity=2093.;

			%ice latent heat of fusion L (J/kg)
			obj.latentheat=3.34*10^5;

			%ice thermal conductivity (W/m/K)
			obj.thermalconductivity=2.4;
			
			%wet ice thermal conductivity (W/m/K)
			obj.temperateiceconductivity=.24;

			%the melting point of ice at 1 atmosphere of pressure in K
			obj.meltingpoint=273.15;

			%rate of change of melting point with pressure (K/Pa)
			obj.beta=9.8*10^-8;

			%mixed layer (ice-water interface) heat capacity (J/kg/K)
			obj.mixed_layer_capacity=3974.;

			%thermal exchange velocity (ice-water interface) (m/s)
			obj.thermal_exchange_velocity=1.00*10^-4;

			%Rheology law: what is the temperature dependence of B with T
			%available: none, paterson and arrhenius
			obj.rheology_law='Paterson';

			% GIA:
			obj.lithosphere_shear_modulus  = 6.7*10^10;  % (Pa)
			obj.lithosphere_density        = 3.32;       % (g/cm^-3)
			obj.mantle_shear_modulus       = 1.45*10^11; % (Pa)
			obj.mantle_density             = 3.34;       % (g/cm^-3)

		end % }}}
		function md = checkconsistency(obj,md,solution,analyses) % {{{
			md = checkfield(md,'fieldname','materials.rho_ice','>',0);
			md = checkfield(md,'fieldname','materials.rho_water','>',0);
			md = checkfield(md,'fieldname','materials.rho_freshwater','>',0);
			md = checkfield(md,'fieldname','materials.mu_water','>',0);
			md = checkfield(md,'fieldname','materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1]);
			md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]);
			md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'Cuffey' 'Paterson' 'Arrhenius' 'LliboutryDuval'});

			if ismember(GiaAnalysisEnum(),analyses),
				md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1);
				md = checkfield(md,'fieldname','materials.lithosphere_density','>',0,'numel',1);
				md = checkfield(md,'fieldname','materials.mantle_shear_modulus','>',0,'numel',1);
				md = checkfield(md,'fieldname','materials.mantle_density','>',0,'numel',1);
			end

		end % }}}
		function disp(obj) % {{{
			disp(sprintf('   Materials:'));

			fielddisplay(obj,'rho_ice','ice density [kg/m^3]');
			fielddisplay(obj,'rho_water','ocean water density [kg/m^3]');
			fielddisplay(obj,'rho_freshwater','fresh water density [kg/m^3]');
			fielddisplay(obj,'mu_water','water viscosity [N s/m^2]');
			fielddisplay(obj,'heatcapacity','heat capacity [J/kg/K]');
			fielddisplay(obj,'thermalconductivity',['ice thermal conductivity [W/m/K]']);
			fielddisplay(obj,'temperateiceconductivity','temperate ice thermal conductivity [W/m/K]');
			fielddisplay(obj,'meltingpoint','melting point of ice at 1atm in K');
			fielddisplay(obj,'latentheat','latent heat of fusion [J/kg]');
			fielddisplay(obj,'beta','rate of change of melting point with pressure [K/Pa]');
			fielddisplay(obj,'mixed_layer_capacity','mixed layer capacity [W/kg/K]');
			fielddisplay(obj,'thermal_exchange_velocity','thermal exchange velocity [m/s]');
			fielddisplay(obj,'rheology_B','flow law parameter [Pa/s^(1/n)]');
			fielddisplay(obj,'rheology_n','Glen''s flow law exponent');
			fielddisplay(obj,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''Cuffey'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']);
			fielddisplay(obj,'lithosphere_shear_modulus','Lithosphere shear modulus [Pa]');
			fielddisplay(obj,'lithosphere_density','Lithosphere density [g/cm^-3]');
			fielddisplay(obj,'mantle_shear_modulus','Mantle shear modulus [Pa]');
			fielddisplay(obj,'mantle_density','Mantle density [g/cm^-3]');
		end % }}}
		function marshall(obj,md,fid) % {{{
			WriteData(fid,'enum',MaterialsEnum(),'data',MaticeEnum(),'format','Integer');
			WriteData(fid,'object',obj,'class','materials','fieldname','rho_ice','format','Double');
			WriteData(fid,'object',obj,'class','materials','fieldname','rho_water','format','Double');
			WriteData(fid,'object',obj,'class','materials','fieldname','rho_freshwater','format','Double');
			WriteData(fid,'object',obj,'class','materials','fieldname','mu_water','format','Double');
			WriteData(fid,'object',obj,'class','materials','fieldname','heatcapacity','format','Double');
			WriteData(fid,'object',obj,'class','materials','fieldname','latentheat','format','Double');
			WriteData(fid,'object',obj,'class','materials','fieldname','thermalconductivity','format','Double');
			WriteData(fid,'object',obj,'class','materials','fieldname','temperateiceconductivity','format','Double');
			WriteData(fid,'object',obj,'class','materials','fieldname','meltingpoint','format','Double');
			WriteData(fid,'object',obj,'class','materials','fieldname','beta','format','Double');
			WriteData(fid,'object',obj,'class','materials','fieldname','mixed_layer_capacity','format','Double');
			WriteData(fid,'object',obj,'class','materials','fieldname','thermal_exchange_velocity','format','Double');
			WriteData(fid,'object',obj,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2);
			WriteData(fid,'data',StringToEnum(obj.rheology_law),'enum',MaterialsRheologyLawEnum(),'format','Integer');

			WriteData(fid,'object',obj,'class','materials','fieldname','lithosphere_shear_modulus','format','Double');
			WriteData(fid,'object',obj,'class','materials','fieldname','lithosphere_density','format','Double','scale',10^3);
			WriteData(fid,'object',obj,'class','materials','fieldname','mantle_shear_modulus','format','Double');
			WriteData(fid,'object',obj,'class','materials','fieldname','mantle_density','format','Double','scale',10^3);
		end % }}}
	end
end