%SNOWPACK class definition
%
%   Usage:
%      snowpack=snowpack();

classdef snowpack
	properties (SetAccess=public) 

		%first, the configuration fields, by category: 
		%snowpack:  %{{{
		snowpack_meas_tss = 0;
		snowpack_enforce_measured_snow_heights = 0;
		snowpack_sw_mode = 0;
		snowpack_incoming_longwave = 0;
		snowpack_height_of_wind_value = 0;
		snowpack_height_of_meteo_values = 0;
		snowpack_neutral = 0;
		snowpack_roughness_length = 0;
		snowpack_number_slopes = 0;
		snowpack_snow_redistribution = 0;
		snowpack_calculation_step_length = 0;
		snowpack_change_bc = 0;
		snowpack_thresh_change_bc = 0;
		snowpack_snp_soil = 0;
		snowpack_soil_flux = 0;
		snowpack_geo_heat = 0;
		snowpack_canopy = 0;
		%}}}
		%snowpackadvanced:  %{{{
		snowpackadvanced_variant = ''; % use 320 kg m-3 for fixed density
		snowpackadvanced_hn_density = '';
		%}}}
		%general:  %{{{
		general_pluginpath = '';
		general_buff_chunk_size = 0;
		general_buff_before = 0;
		%}}}
		%input {{{
		input_coordsys = '';
		input_coordparam = '';
		input_time_zone = 0;
		input_meteo = '';
		input_meteopath = '';
		input_station1 = '';
		input_snowfile1 = '';
		%}}}
		%output {{{
		output_coordsys = '';
		output_coordparam = '';
		output_time_zone = 0;
		output_meteopath = '';
		output_experiment = '';
		output_ts_write = 0;
		output_ts_start = 0;
		output_ts_days_between = 0;
		output_profile = '';
		output_prof_write = 0;
		output_prof_start = 0;
		output_prof_days_between = 0;
		%}}}
		%interpolations1d %{{{
		interpolations1d_window_size = 0; %that is 5 d and 2 h; 1 d = 86400
		interpolations1d_hnw_resample = '';
		interpolations1d_hs_resample = '';
		interpolations1d_tsg_resample = '';
		interpolations1d_rho_hn_resample = '';
		interpolations1d_vw_resample = '';
		interpolations1d_vw_args = '';
		%}}}
		%filters {{{
		filters_ta_filter1 = '';
		filters_ta_arg1 = NaN;
		filters_rh_filter1 = '';
		filters_rh_arg1 = NaN;
		filters_rh_filter2 = '';
		filters_rh_arg2 = NaN;
		filters_iswr_filter1 = '';
		filters_iswr_arg1 = NaN;
		filters_iswr_filter2 = '';
		filters_iswr_arg2 = NaN;
		filters_rswr_filter1 = '';
		filters_rswr_arg1 = NaN;
		filters_rswr_filter2 = '';
		filters_rswr_arg2 = NaN;

		%for ta between 190 and 280 k;
		filters_ilwr_filter1 = '';
		filters_ilwr_arg1 = NaN;
		filters_ilwr_filter2 = '';
		filters_ilwr_arg2 = NaN;
		filters_tss_filter1 = '';
		filters_tss_arg1 = NaN;
		filters_tsg_filter1 = '';
		filters_tsg_arg1 = NaN;
		filters_vw_filter1 = '';
		filters_vw_arg1 = NaN;
		filters_vw_filter2 = '';
		filters_vw_arg2 = NaN;
		%}}}

	end
	methods
		function obj = snowpack(varargin) % {{{
			switch nargin
				case 0
					obj=setdefaultparameters(obj);
				case 1
					inputstruct=varargin{1};
					list1 = properties('snowpack');
					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) % {{{
	
		%snowpack:  %{{{
		obj.snowpack_meas_tss = 1;
		obj.snowpack_enforce_measured_snow_heights = 0;
		obj.snowpack_sw_mode = 0;
		obj.snowpack_incoming_longwave = 1;
		obj.snowpack_height_of_wind_value = 12.;
		obj.snowpack_height_of_meteo_values = 12.;
		obj.snowpack_neutral = 0;
		obj.snowpack_roughness_length = 0.002;
		obj.snowpack_number_slopes = 1;
		obj.snowpack_snow_redistribution = 1;
		obj.snowpack_calculation_step_length = 15.0;
		obj.snowpack_change_bc = 0;
		obj.snowpack_thresh_change_bc = -1.0;
		obj.snowpack_snp_soil = 0;
		obj.snowpack_soil_flux = 0;
		obj.snowpack_geo_heat = 0.06;
		obj.snowpack_canopy = 0;
		%}}}
		%snowpackadvanced:  %{{{
		obj.snowpackadvanced_variant = 'ANTARCTICA'; % use 320 kg m-3 for fixed density
		obj.snowpackadvanced_hn_density = 'EVENT';
		%}}}
		%general:  %{{{
		obj.general_pluginpath = '/usr/local/lib/meteoio/plugins/';
		obj.general_buff_chunk_size = 90;
		obj.general_buff_before = 1.5;
		%}}}
		%input {{{
		obj.input_coordsys = 'ch1903';
		obj.input_coordparam = 'null';
		obj.input_time_zone = 8;
		obj.input_meteo = 'smet';
		obj.input_meteopath = './input';
		obj.input_station1 = 'domec.smet';
		obj.input_snowfile1 = 'domec.sno';
		%}}}
		%output {{{
		obj.output_coordsys = 'ch1903';
		obj.output_coordparam = 'null';
		obj.output_time_zone = 8;
		obj.output_meteopath = './output';
		obj.output_experiment = 'smet';
		obj.output_ts_write = 1;
		obj.output_ts_start = 0.0;
		obj.output_ts_days_between = 0.04166667;
		obj.output_profile = 'ascii';
		obj.output_prof_write = 1;
		obj.output_prof_start = 0.0;
		obj.output_prof_days_between = 0.04166667;
		%}}}
		%interpolations1d %{{{
		obj.interpolations1d_window_size = 439200 %that is 5 d and 2 h; 1 d = 86400
		obj.interpolations1d_hnw_resample = 'none';
		obj.interpolations1d_hs_resample = 'linear';
		obj.interpolations1d_tsg_resample = 'linear';
		obj.interpolations1d_rho_hn_resample = 'none';
		obj.interpolations1d_vw_resample = 'nearest_neighbour';
		obj.interpolations1d_vw_args = 'extrapolate';
		%}}}
		%filters {{{
		obj.filters_ta_filter1 = 'min_max';
		obj.filters_ta_arg1 = [190 280];
		obj.filters_rh_filter1 = 'min_max';
		obj.filters_rh_arg1 = [0.01 1.2];
		obj.filters_rh_filter2 = 'min_max';
		obj.filters_rh_arg2 = {'soft' 0.01 1.0};
		obj.filters_iswr_filter1 = 'min_max';
		obj.filters_iswr_arg1 = [-10 1500];
		obj.filters_iswr_filter2 = 'min_max';
		obj.filters_iswr_arg2 = {'soft' 0 1500};
		obj.filters_rswr_filter1 = 'min_max';
		obj.filters_rswr_arg1 = [-10 1500];
		obj.filters_rswr_filter2 = 'min_max';
		obj.filters_rswr_arg2 = {'soft' 0 1500};

		%for ta between 190 and 280 k;
		obj.filters_ilwr_filter1 = 'min_max';
		obj.filters_ilwr_arg1 = [30 355];
		obj.filters_ilwr_filter2 = 'min_max';
		obj.filters_ilwr_arg2 = {'soft' 35 350};
		obj.filters_tss_filter1 = 'min_max';
		obj.filters_tss_arg1 = [180 275];
		obj.filters_tsg_filter1 = 'min_max';
		obj.filters_tsg_arg1 = [200 275];
		obj.filters_vw_filter1 = 'min_max';
		obj.filters_vw_arg1 = [-2 70];
		obj.filters_vw_filter2 = 'min_max';
		obj.filters_vw_arg2 = {'soft' 0 50};
		%}}}

		end % }}}
		function md = checkconsistency(obj,md,solution,analyses) % {{{
			%snowpack:  %{{{
			md=checkfield(md,'snowpack.snowpack_meas_tss','values',[0 1]);
			md=checkfield(md,'snowpack.snowpack_enforce_measured_snow_heights','values',[0 1]);
			md=checkfield(md,'snowpack.snowpack_sw_mode','values',[0 1 2]);
			md=checkfield(md,'snowpack.snowpack_incoming_longwave','values',[0 1]);
			md=checkfield(md,'snowpack.snowpack_height_of_wind_value,'>=',0);
			md=checkfield(md,'snowpack.snowpack_height_of_meteo_values','>=',0);
			md=checkfield(md,'snowpack.snowpack_neutral','values',[-1 0 1]);
			md=checkfield(md,'snowpack.snowpack_roughness_length','>=',0);
			md=checkfield(md,'snowpack.snowpack_number_slopes','values',[1 3 5 9]);
			md=checkfield(md,'snowpack.snowpack_snow_redistribution','values',[0 1]);
			md=checkfield(md,'snowpack.snowpack_calculation_step_length','>',0);
			md=checkfield(md,'snowpack.snowpack_change_bc','values',[0 1]);
			md=checkfield(md,'snowpack.snowpack_thresh_change_bc','<=',0);
			md=checkfield(md,'snowpack.snowpack_snp_soil','values',[0 1]);
			md=checkfield(md,'snowpack.snowpack_soil_flux,'values',[0 1]);
			md=checkfield(md,'snowpack.snowpack_geo_heat,'>=',0);
			md=checkfield(md,'snowpack.snowpack_canopy,'values',[0 1]);
			%}}}
			%snowpackadvanced:  %{{{
			md=checkfield(md,'snowpack.snowpackadvanced_variant = 'ANTARCTICA'; % use 320 kg m-3 for fixed density
			md=checkfield(md,'snowpack.snowpackadvanced_hn_density = 'EVENT';
			%}}}
			%general:  %{{{
			md=checkfield(md,'snowpack.general_pluginpath = '/usr/local/lib/meteoio/plugins/';
			md=checkfield(md,'snowpack.general_buff_chunk_size = 90;
			md=checkfield(md,'snowpack.general_buff_before = 1.5;
			%}}}
			%input {{{
			md=checkfield(md,'snowpack.input_coordsys = 'ch1903';
			md=checkfield(md,'snowpack.input_coordparam = 'null';
			md=checkfield(md,'snowpack.input_time_zone = 8;
			md=checkfield(md,'snowpack.input_meteo = 'smet';
			md=checkfield(md,'snowpack.input_meteopath = './input';
			md=checkfield(md,'snowpack.input_station1 = 'domec.smet';
			md=checkfield(md,'snowpack.input_snowfile1 = 'domec.sno';
			%}}}
			%output {{{
			md=checkfield(md,'snowpack.output_coordsys = 'ch1903';
			md=checkfield(md,'snowpack.output_coordparam = 'null';
			md=checkfield(md,'snowpack.output_time_zone = 8;
			md=checkfield(md,'snowpack.output_meteopath = './output';
			md=checkfield(md,'snowpack.output_experiment = 'smet';
			md=checkfield(md,'snowpack.output_ts_write = 1;
			md=checkfield(md,'snowpack.output_ts_start = 0.0;
			md=checkfield(md,'snowpack.output_ts_days_between = 0.04166667;
			md=checkfield(md,'snowpack.output_profile = 'ascii';
			md=checkfield(md,'snowpack.output_prof_write = 1;
			md=checkfield(md,'snowpack.output_prof_start = 0.0;
			md=checkfield(md,'snowpack.output_prof_days_between = 0.04166667;
			%}}}
			%interpolations1d %{{{
			md=checkfield(md,'snowpack.interpolations1d_window_size = 439200 %that is 5 d and 2 h; 1 d = 86400
			md=checkfield(md,'snowpack.interpolations1d_hnw_resample = 'none';
			md=checkfield(md,'snowpack.interpolations1d_hs_resample = 'linear';
			md=checkfield(md,'snowpack.interpolations1d_tsg_resample = 'linear';
			md=checkfield(md,'snowpack.interpolations1d_rho_hn_resample = 'none';
			md=checkfield(md,'snowpack.interpolations1d_vw_resample = 'nearest_neighbour';
			md=checkfield(md,'snowpack.interpolations1d_vw_args = 'extrapolate';
			%}}}
			%filters {{{
			md=checkfield(md,'snowpack.filters_ta_filter1 = 'min_max';
			md=checkfield(md,'snowpack.filters_ta_arg1 = [190 280];
			md=checkfield(md,'snowpack.filters_rh_filter1 = 'min_max';
			md=checkfield(md,'snowpack.filters_rh_arg1 = [0.01 1.2];
			md=checkfield(md,'snowpack.filters_rh_filter2 = 'min_max';
			md=checkfield(md,'snowpack.filters_rh_arg2 = {'soft' 0.01 1.0};
			md=checkfield(md,'snowpack.filters_iswr_filter1 = 'min_max';
			md=checkfield(md,'snowpack.filters_iswr_arg1 = [-10 1500];
			md=checkfield(md,'snowpack.filters_iswr_filter2 = 'min_max';
			md=checkfield(md,'snowpack.filters_iswr_arg2 = {'soft' 0 1500};
			md=checkfield(md,'snowpack.filters_rswr_filter1 = 'min_max';
			md=checkfield(md,'snowpack.filters_rswr_arg1 = [-10 1500];
			md=checkfield(md,'snowpack.filters_rswr_filter2 = 'min_max';
			md=checkfield(md,'snowpack.filters_rswr_arg2 = {'soft' 0 1500};

			%for ta between 190 and 280 k;
			md=checkfield(md,'snowpack.filters_ilwr_filter1 = 'min_max';
			md=checkfield(md,'snowpack.filters_ilwr_arg1 = [30 355];
			md=checkfield(md,'snowpack.filters_ilwr_filter2 = 'min_max';
			md=checkfield(md,'snowpack.filters_ilwr_arg2 = {'soft' 35 350};
			md=checkfield(md,'snowpack.filters_tss_filter1 = 'min_max';
			md=checkfield(md,'snowpack.filters_tss_arg1 = [180 275];
			md=checkfield(md,'snowpack.filters_tsg_filter1 = 'min_max';
			md=checkfield(md,'snowpack.filters_tsg_arg1 = [200 275];
			md=checkfield(md,'snowpack.filters_vw_filter1 = 'min_max';
			md=checkfield(md,'snowpack.filters_vw_arg1 = [-2 70];
			md=checkfield(md,'snowpack.filters_vw_filter2 = 'min_max';
			md=checkfield(md,'snowpack.filters_vw_arg2 = {'soft' 0 50};
			%}}}

		end % }}}
		function disp(obj) % {{{

			disp(sprintf('   StressBalance solution parameters:'));

			disp(sprintf('\n      %s','Convergence criteria:'));
			fielddisplay(obj,'restol','mechanical equilibrium residual convergence criterion');
			fielddisplay(obj,'reltol','velocity relative convergence criterion, NaN: not applied');
			fielddisplay(obj,'abstol','velocity absolute convergence criterion, NaN: not applied');
			fielddisplay(obj,'isnewton','0: Picard''s fixed point, 1: Newton''s method, 2: hybrid');
			fielddisplay(obj,'maxiter','maximum number of nonlinear iterations');
			fielddisplay(obj,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)');

			disp(sprintf('\n      %s','boundary conditions:'));
			fielddisplay(obj,'spcvx','x-axis velocity constraint (NaN means no constraint) [m/yr]');
			fielddisplay(obj,'spcvy','y-axis velocity constraint (NaN means no constraint) [m/yr]');
			fielddisplay(obj,'spcvz','z-axis velocity constraint (NaN means no constraint) [m/yr]');

			disp(sprintf('\n      %s','Rift options:'));
			fielddisplay(obj,'rift_penalty_threshold','threshold for instability of mechanical constraints');
			fielddisplay(obj,'rift_penalty_lock','number of iterations before rift penalties are locked');

			disp(sprintf('\n      %s','Penalty options:'));
			fielddisplay(obj,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset');
			fielddisplay(obj,'vertex_pairing','pairs of vertices that are penalized');

			disp(sprintf('\n      %s','Other:'));
			fielddisplay(obj,'shelf_dampening','use dampening for floating ice ? Only for FS model');
			fielddisplay(obj,'FSreconditioning','multiplier for incompressibility equation. Only for FS model');
			fielddisplay(obj,'referential','local referential');
			fielddisplay(obj,'loadingforce','loading force applied on each point [N/m^3]');
			fielddisplay(obj,'requested_outputs','additional outputs requested');

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

			yts=365.0*24.0*3600.0;

			WriteData(fid,'object',obj,'class','snowpack','fieldname','spcvx','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
			WriteData(fid,'object',obj,'class','snowpack','fieldname','spcvy','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
			WriteData(fid,'object',obj,'class','snowpack','fieldname','spcvz','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1);
			WriteData(fid,'object',obj,'class','snowpack','fieldname','restol','format','Double');
			WriteData(fid,'object',obj,'class','snowpack','fieldname','reltol','format','Double');
			WriteData(fid,'object',obj,'class','snowpack','fieldname','abstol','format','Double');
			WriteData(fid,'object',obj,'class','snowpack','fieldname','isnewton','format','Integer');
			WriteData(fid,'object',obj,'class','snowpack','fieldname','FSreconditioning','format','Double');
			WriteData(fid,'object',obj,'class','snowpack','fieldname','viscosity_overshoot','format','Double');
			WriteData(fid,'object',obj,'class','snowpack','fieldname','maxiter','format','Integer');
			WriteData(fid,'object',obj,'class','snowpack','fieldname','shelf_dampening','format','Integer');
			WriteData(fid,'object',obj,'class','snowpack','fieldname','vertex_pairing','format','DoubleMat','mattype',3);
			WriteData(fid,'object',obj,'class','snowpack','fieldname','penalty_factor','format','Double');
			WriteData(fid,'object',obj,'class','snowpack','fieldname','rift_penalty_lock','format','Integer');
			WriteData(fid,'object',obj,'class','snowpack','fieldname','rift_penalty_threshold','format','Integer');
			WriteData(fid,'object',obj,'class','snowpack','fieldname','referential','format','DoubleMat','mattype',1);
			WriteData(fid,'object',obj,'class','snowpack','fieldname','requested_outputs','format','DoubleMat','mattype',3);
			WriteData(fid,'data',obj.loadingforce(:,1),'format','DoubleMat','mattype',1,'enum',LoadingforceXEnum);
			WriteData(fid,'data',obj.loadingforce(:,2),'format','DoubleMat','mattype',1,'enum',LoadingforceYEnum);
			WriteData(fid,'data',obj.loadingforce(:,3),'format','DoubleMat','mattype',1,'enum',LoadingforceZEnum);
		end % }}}
	end
end