function [rt_g m]=cielothermal_core(m,params,thermalparams,inputs,analysis)
%INPUT function rt_g =cielothermal_core(m,params,thermalparams,inputs,analysis)

	%check inputs
	if isempty(inputs),
		clear inputs;
	end
	
	%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

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

	dt=inputs.dt;

	disp(sprintf('%s','   starting direct shooting method'));

	%Recover pressure from inputs
	pressure=inputs.pressure;

	%Optional inputs
	inputs.beta=thermalparams.beta; 
	inputs.meltingpoint=thermalparams.meltingpoint;
	inputs.mintemp=thermalparams.mintemp;

	%intialize lst
	lst=m.lst;
	
	count=1;

	while(~converged),

		%system matrices
		[rK_gg , rp_g , rdK_gg]=Emg(m.bgpdt,m.bgpdtb, m.est,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 ); 
	
		%Cleanup memory
		rK_gg=IMdb('drop rK_gg');
		if ~isempty(rdK_gg), rdK_gg=IMdb('drop 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 );

		%Cleanup memory
		[rp_g rK_fs]=IMdb('drop rp_g rK_fs');

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

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

		%Deal with penalty loads
		[new_lst,converged,num_unstable_constraints]=CieloPenaltyConstraints(m.bgpdt,m.bgpdtb,m.est,lst, struct('pressure',pressure, 'dt',dt, 'temperature',rt_g,'beta',thermalparams.beta,'melting_point',thermalparams.meltingpoint),analysis);

		if ~converged,
			disp(sprintf('   %s %i','#unstable constraints ',num_unstable_constraints));
			if num_unstable_constraints<thermalparams.min_thermal_constraints,
				converged=1;
			end
			if count==2,
				disp('Carefull! Screwed up solutoin on purpose!');
				converged=1;
			end
		end

		%Set lst to new_lst, but do not wipeout m.lst at the same time! (we can't modify the input m.lst)
		lst=Equiv(lst,new_lst);

		count=count+1;
	end

	%don't forget to output loads
	m.lst=lst;
	
end
