%Dynamic heat transfer
	clear all
	close all

%different indentation level for comments for automated text folding
%start of simulation 
	runorjustplot=input('Run or just plot? (1 for run, 2 for just plot) ');
	if runorjustplot==1
		simulationset=input(['simulationset=? (1 for 15 ' char(176) 'C/min heating, 15-mm-thick specimen; 2 for 2 ' char(176) 'C/min, 1-mm-thick specimen. 2 is the situation reflecting the experiment conditions)']);
		initialtempexperiment=25; %initial temperature of the problem (experiment) to be solved
		initialtempsimulate=initialtempexperiment; %always needs to be equal to initialtempexperiment for correct simulation. initialtempsimulate larger than initialtempexperiment only for testing.
		endtemp=200; %temperature at the end of the simulation
		speedupfactor=1; %needs to be equal to 1 for correct simulation, other values only for testing
		if simulationset==1	
			heatingrate=15/60; %degrees per second
		else
			heatingrate=2/60; %degrees per second
		end
		heatingrate=heatingrate*speedupfactor; %heating rate after altered with speedupfactor, speedupfactor=1 for correct simulation, other values for testing purposes only.
		totaltime=(endtemp-initialtempsimulate)/heatingrate; %s
		if simulationset==1
			totaldepth=15/1000; %m (15 mm or 15/1000 m)
		else
			totaldepth=1/1000; %m (1 mm or 1/1000 m)
		end
		ndepthpoint=101; %number of nodes
		w=0.01; %specimen width in meters
		k=@(T)0.1820385*(1+0.002*(T-initialtempexperiment)); %thermal conductivity (function of temperature); %https://doi.org/10.1016/0960-8524(91)90205-X


		c=@(T)(1114+4.86*T);	%specific heat capacity (function of temperature); Kollmann, F., Coˆte´, W.A., Jr. (1968) Principles of Wood Science and Technology.

		rho=425; %density; kg/m3
		alpha_transfer=@(T)k(T)/(c(T)*rho); %alpha for heat diffusivity (function of temperature)
		alpha_expand=34.6*10^-6; %alpha for thermal expansion; Weatherwax, R. C., & Stamm, A. J. (1956). The coefficients of thermal expansion of wood and wood products.

		dx=totaldepth/(ndepthpoint-1); %distance between nodes
		dt=dx^2/(2*max(alpha_transfer(25),alpha_transfer(200))); %time step size
		ntimepoint=ceil(totaltime/dt+1); %number of discretised time steps
		u=zeros(ndepthpoint+1,2); %temperature (only keeping two time steps to save memory). we keep track of temperature every several timesteps just for plotting
		epsilonT=zeros(ndepthpoint,1); %thermal expansion strain
		Ftarget=0.9; %constantly-applied tension; N
		%Ftarget=0;
		Einitial=820*1000000; %Elastic modulus at the start of experiment; Pa
		epsilon_initial=Ftarget/w/totaldepth/Einitial;
		E=zeros(1,ndepthpoint); %Preallocating elastic moduli values at nodes
		EtimesepsilonT=zeros(1,ndepthpoint);
		dF=zeros(1,ndepthpoint);
		depth=zeros(1,ndepthpoint);
		airtemp=zeros(1,ntimepoint+1);
		deltaT=zeros(1,ntimepoint+1);
		time=zeros(1,ntimepoint+1);
		for idepthpoint=1:ndepthpoint
			 depth(idepthpoint)=(idepthpoint-1)*dx;
		end
		for itimepoint=1:ntimepoint+1
			 time(itimepoint)=(itimepoint-1)*dt;
			 airtemp(itimepoint)=initialtempsimulate+(itimepoint-1)/(ntimepoint-1)*(endtemp-initialtempsimulate);
		end


		%start of plot variable declarations
		ntimepointplot=11; %number of time points grabbed for surface plotting in 3D
		ntimepointplot2=201; %number of time points grabbed for line plotting in 2D
		epsilont=zeros(1,ntimepointplot2-1); %total strain, preallocated for grabbing values for line plot
		deltaTplot=zeros(1,ntimepointplot2-1); %difference between highest and lowest temperature at a given time, preallocated for grabbing values for line plot
		ndepthpointplot=11; %number of depth points for the surface plot in 3D
		dtplot=totaltime/(ntimepointplot-1);
		dtplot2=totaltime/(ntimepointplot2-1);
		dxplot=totaldepth/(ndepthpointplot-1);
		timepointplot=zeros(ndepthpointplot-1,ntimepointplot-1);
		timepointplot2=zeros(1,ntimepointplot2);
		depthpointplot=zeros(ndepthpointplot-1,ntimepointplot-1);
		timepointplotregression=zeros(1,ntimepointplot);
		depthpointplotregression=zeros(1,ndepthpointplot);
		upointplot=zeros(ndepthpointplot-1,ntimepointplot-1);
		upointplot2=zeros(ndepthpoint,ntimepointplot2-1);
		airtempplot=zeros(1,ntimepointplot2-1);
		epsilontplot=zeros(1,ntimepointplot2-1);
		%end of plot variable declarations

		%start if setting boundaries
		for idepthpoint=2:ndepthpoint-1
			u(idepthpoint,1)=initialtempsimulate;
		end
		%end of setting boundaries

		%start of explicit time integration
		itimepointplot=1;
		itimepointplot2=1;
		u(1,1)=airtemp(1);
		u(ndepthpoint,1)=airtemp(1);
		for itimepoint=1:ntimepoint
			u(1,2)=airtemp(itimepoint+1);
			u(ndepthpoint,2)=airtemp(itimepoint+1);
			for idepthpoint=2:ndepthpoint-1
				d2uperdx2=(u(idepthpoint+1,1)+u(idepthpoint-1,1)-(2*u(idepthpoint,1)))/dx^2;
				u(idepthpoint,2)=u(idepthpoint,1)+alpha_transfer(u(idepthpoint,1))*d2uperdx2*dt;
			end
			u(:,1)=u(:,2);
			if (itimepointplot-1)*dtplot<=(itimepoint-1)*dt
				idepthpointplot=1;
				for idepthpoint=1:ndepthpoint
					if (idepthpointplot-1)*dxplot<=(idepthpoint-1)*dx
						timepointplot(idepthpointplot,itimepointplot)=(itimepoint-1)*dt;
						depthpointplot(idepthpointplot,itimepointplot)=(idepthpoint-1)*dx;
						upointplot(idepthpointplot,itimepointplot)=u(idepthpoint,1);
						idepthpointplot=idepthpointplot+1;
					end
				end
				itimepointplot=itimepointplot+1;
				disp([num2str(itimepoint/ntimepoint*100) '% complete']);
			end
			if (itimepointplot2-1)*dtplot2<=(itimepoint-1)*dt
				airtempplot(itimepointplot2)=airtemp(itimepoint);
				for idepthpoint=1:ndepthpoint
					epsilonT(idepthpoint)=alpha_expand*(u(idepthpoint,1)-initialtempexperiment);
					E(idepthpoint)=Einitial*(1-(u(idepthpoint,1)-initialtempexperiment)/(300-initialtempexperiment));
					EtimesepsilonT(idepthpoint)=E(idepthpoint)*epsilonT(idepthpoint);
				end
				epsilont(itimepointplot2)=(Ftarget/w+trapz(depth,EtimesepsilonT))/(trapz(depth,E));
				deltaTplot(itimepointplot2)=u(1,1)-u(round((ndepthpoint-1)/2),1);
				itimepointplot2=itimepointplot2+1;
			end
		end
		%end of explicit time integration
	else
		disp('Please load the run data')
		uiload;
		runorjustplot=2; %since the loaded .mat file will have runorjustplot=1, this is to correct it
	end
%end of simulation

%start of outputs
	alpha_expand2=(epsilont(itimepointplot2-1)-epsilon_initial)/(airtempplot(itimepointplot2-1)-initialtempexperiment);
	alpha_expand_error=abs(alpha_expand2-alpha_expand)/alpha_expand;
	disp(['alpha_expand_error=' num2str(alpha_expand_error)]);
	
	close all;
	ifigure=0;

	%start of 3d plot arrangement
	ifigure=ifigure+1;
	figure(ifigure);
	surf(depthpointplot*1000,timepointplot/60,upointplot);
	xlabel(['Specimen Depth (mm)']);
	ylabel(['Time (minutes)']);
	zlabel(['Temperature (' char(176) 'C)']);
	colormap jet;
	title(['Temperature Development' char(10) 'Throughout Specimen Thickness']);
	pause(0.001);
	set(figure(ifigure),'Position',[0 0 400 300]);
	%end of 3d plot arrangement

	%start of 1d epsilont vs temp plot
	ifigure=ifigure+1;
	figure(ifigure);
	epsilonapparent=epsilont-epsilon_initial;
	epsilontheoretical=alpha_expand*(airtempplot(itimepointplot2-1)-initialtempexperiment);
	strainerror=(epsilonapparent(itimepointplot2-1)-epsilontheoretical)/epsilontheoretical; %should be the same as alpha_expand_error
	disp(['epsilonapparent(t_final)=' num2str(epsilonapparent(itimepointplot2-1))]);
	disp(['epsilontheoretical=' num2str(epsilontheoretical)]);
	disp(['strainerror=' num2str(strainerror)]);
	plot(airtempplot,epsilonapparent);
	title('Apparent Swell Strain');
	xlabel(['Oven Temperature (' char(176) 'C)']);
	ylabel('Apparent Strain (mm/mm)');
	set(figure(ifigure),'Position',[0 0 400 300]);
	%end of 1d epsilont vs time plot

	%start of 1d deltaT vs temp plot
	ifigure=ifigure+1;
	figure(ifigure);
	plot(airtempplot,deltaTplot);
	title('Highest - Lowest Temperature');
	xlabel(['Oven Temperature (' char(176) 'C)']);
	ylabel(['Max Temp. Difference (' char(176) 'C )']);
	set(figure(ifigure),'Position',[0 0 400 300]);
	%end of 1d deltaT vs temp plot
%end of outputs

%start of saves
	if runorjustplot==1
		save results.mat; %saves the results as resuts.mat in case need to plot again later. change the .mat file name into something else to have differen save file. Remember to also change the loaded .mat file at the start of the script if runorjustplot is set to 2
	end
%end of saves