function par = makeDefaultParameters(rainType, baseCatchment)
	% rainType: "pulse", "poisson", or "observed"
	% baseCatchment: specify a catchment name to use its empirical parameters

	% define catchment combination if not specified
	if ~exist("baseCatchment", "var")
		baseCatchment = "HubbardBrook25";
	end

	% define parameters
	par = defineDefaults(rainType);
	par = updateByCatchment(par, baseCatchment);
	par = makeDependentParameters(par);
end

function par = defineDefaults(rainType)
	% network length parameters
	par.domainL.A = 5;		 % [km2]  catchment area
	par.domainL.phiC = 0.75; % [-]	  non-perennial fraction
	par.domainL.Dg = 3.68;	 % [km-1] geomorphic drainage density
	par.domainL.dL = 50;	 % [m]	  resolution of length domain for solving f(tau, L)
	% par.domainL.Lmax;		 % [m]	  maximum network length
	%par.domainL.Lmin = 8000;		 % [m]	  minimum network length
	% par.domainL.Lrange;	 % [m]	  range of variability of network length

	% alpha function
	par.alpha.beta = 160; %120;
	% par.alpha.fun = @(L, tau) par.alpha.beta * tau ./ (par.alpha.beta * tau + L);

	% boundary condition for f
	par.f0.k = 0.018; %1/24; % [1/h] inverse of mean response time for exponential iuh as boundary cond of f
	par.f0.theta = 0.803; % [-] shape of gamma IUH at Lmax (theta=1 -> exponential); use with analyticalSolution
	par.f0.uh = 20.2; %1e-02 * 3600; % [m/h] water celerity in hyporheic zone of dry network;
	% par.f0.L0; % [m] where the initial condition is applied
	% par.f0.fun = @(tau) par.f0.k * exp(-par.f0.k * tau); % [1/h] boundary condition (iuh at L = L0)

	% length-discharge relation. Default is hubbard brook
	% see SI of this paper for empirical values: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2024WR038732
	par.LQ.k = 0.767; % [-] shape parameter of gamma function
	par.LQ.theta = 4.45/24; % [mm/h] scale parameter of the distribution
	% par.LQ.fun = @(Q) par.domainL.Lmin + par.domainL.Lrange * gammainc(Q / par.LQ.theta, par.LQ.k);

	% time parameters (default; case "observed" may override)
	par.domainT.D = 1095*24; % [h] total duration of simulation
	par.domainT.dt = 1/8;    % [h] fixed time step
	% par.domainT.Nt; %      [-] number of time steps

	% tau parameters
	par.domainTau.tauMax = 365*24; % [h] make sure it's long enough such that IUH is always complete
	% par.domainTau.dtau;   [h] resolution in tau = resolution in t to simplify convolution
	% par.domainTau.Ntau;   [-] number of steps in tau

	% rainfall parameters
	switch lower(rainType)
		case "pulse" % define rainfall (single constant event)
			par.rain.type = @generateRainfallPulse;
			par.rain.start = 24; % [h] starting hour of rainfall
			par.rain.duration = 1; % [h]duration of rainfall
			par.rain.volume = 40; % [mm] total rain volume (100/24*rain.duration = 100 mm/d)
			par.rain.shape = "rectangular"; % "rectangular", "triangular"
		case "poisson"
			% define rainfall (poisson)
			par.rain.type = @generateRainfallPoisson;
			par.rain.alpha = 15; % [mm/d] avg daily rain depth (exponentially distributed)
			par.rain.lambda = 0.3; % [1/d] avg daily rain frequency (i.e., probability of rainfall for any given day)
			par.rain.downscale = "conservative"; % downscaling type: "linear", "triangular", "conservative"
			par.rain.downscaleLimit = 80/24; % [mm/h] intensity limit after downscaling
			par.rain.downscaleBeta = 1; % beta distrib param to extract downscaling factors. 1 = uniform distribution; >1 = peaked around middle
		case "observed"
			par.rain.type = @generateObservedRainfall;
			par.rain.fileInput = 'meteoData.xlsx';
			par.rain.spdFile = 144; % [1/d] passi nel file (144 = ogni 10 min)

			par.rain.et.latitude = 45.85;    % [°N] Valfredda
			par.rain.et.elevation = 2100;  % [m] per Penman-Monteith (futuro)
			par.rain.et.Kc = 0.83;
			par.rain.et.ETsw_frac = 0.05;

			% Suolo: valori in cm/d e cm (convertiti in mm/h in generateObservedRainfall)
			par.rain.soil.Ks = 201.7;  % 873.6;  
			par.rain.soil.b = 1.6; %5.04;
			par.rain.soil.n = 0.433;
			par.rain.soil.Zr = 30;
			par.rain.soil.s0 = 0.484; %0.50;
			par.rain.soil.sw = 0.15; %0.10;
			par.rain.soil.sstar = 0.55; %0.47

			par.domainT.dt = 1;    % [h] risoluzione meteo e bilancio suolo
			par.domainT.D = 68*24; % [h] ~3 mesi meteo + maggio intero (67.5*24 tagliava le ultime ore)

			par.obs.fileInput = 'measuredDischarge.xlsx';
			par.obs.spdFile = 144;   % passi/giorno nel file (10 min)
			par.obs.units = 'm3s';   % unità nel file
			par.obs.plotLastDays = 30; % giorni da plottare vs Q osservata (fine serie)

		otherwise
			error("rainType must be 'pulse', 'poisson', or 'observed'");
	end
	
	par.rain.baseflow = 0.2; % [mm/h] constant baseflow to be added to flowrate 0.1076
end

function par = updateByCatchment(par, baseCatchment)
	% allocate defaults
	catchments = table( ...
		'Size', [0 6], ...
		'VariableTypes', ["string", repmat("doublenan", 1, 5)], ...
		'VariableNames', ["name",				"A",	"Dg",	"phiC", "LQ_k", "LQ_theta"] ...
	);

	% set up defaults table
	catchments(end+1, :) = {"HubbardBrook25",	0.25,	8.56,	0.79,	0.767,	4.43/24};
	catchments(end+1, :) = {"Valfredda",		5.21,	3.68,	0.75,	0.175,	61.9/24}; 
	catchments(end+1, :) = {"PovertyCreek",		0.35,	5.13,	0.99,	1.43,	19.2/24};
	catchments(end+1, :) = {"DevonA",			0.35,	6.91,	0.84,	0.321,	79.3/24};
%velfredda prima 0.175, 61.9/24
	% find catchment
	id = find(catchments.name == baseCatchment, 1);

	% check specified one exists
	if isempty(id)
		id = 1;
		warning("Specified catchment is not defined. Falling back to default parameters (%s)", catchments.name(id));
	end

	% update params
	par.baseCatchment = catchments.name(id);
	par.domainL.A = catchments.A(id);
	par.domainL.Dg = catchments.Dg(id);
	par.domainL.phiC = catchments.phiC(id);
	par.LQ.k = catchments.LQ_k(id);
	par.LQ.theta = catchments.LQ_theta(id);
end