function rmelting_g =cielomelting_core_linear(m,params,thermalparams,inputs,analysis)
%INPUT function rmelting_g =cielomelting_core_linear(m,params,thermalparams,inputs,analysis)

	%recover flag_y_s from m.y_s
	if m.y_s.header.M==0,
		flag_y_s= 0;
	else
		flag_y_s= 1;
	end

	%stiffness and load generation only:
	params.kflag=1; params.pflag=1;
	params.ktflag=0;

	%optional inputs
	inputs.beta=thermalparams.beta;
	inputs.meltingpoint=thermalparams.meltingpoint;
	inputs.latentheat=thermalparams.latentheat;
	inputs.heatcapacity=thermalparams.heatcapacity;
	inputs.penalty=thermalparams.penalty;

	%system matrices
	[rK_gg , rp_g , rdK_gg]=Emg(m.bgpdt,m.bgpdtb, m.est,m.lst,m.ept,m.mpt,m.geom3,params,inputs,analysis);

	%Reduce tangent matrix from g size to f size
	[rK_ff, rK_fs] = Reducematrixfromgtof( rK_gg, m.G_mn, flag_y_s ); 

	if params.lowmem, 
		[rK_gg,rdK_gg]=IMdb('drop rK_gg rdK_gg');
	end
	
	%Reduce load from g size to f size
	[rp_f] = Reducerightside( rp_g, m.G_mn, rK_fs, m.y_s, flag_y_s );
	if params.lowmem,
		[rp_g rK_fs]=IMdb('drop rp_g rK_fs');
	end

	%Solve	
	[rt_f]=Solver(rK_ff,rp_f,{},params.solverstring);

	if params.lowmem,
		[rK_ff rp_f]=IMdb('drop rK_ff rp_f');
	end
	
	%Merge back to g set
	rmelting_g= Mergesolvec( rt_f, m.G_mn, m.y_s ); 
	
	if (params.lowmem && count>2),
		rt_f=IMdb('drop rt_f ');
	end
