From 918907c2e89dc8b46163f16cbbd56946774cadd5 Mon Sep 17 00:00:00 2001
From: Robert Chen <fernjager@gmail.com>
Date: Sun, 1 Dec 2013 22:51:46 -0800
Subject: [PATCH] Matlab simulation file

---
 sim/TurkeySimulation.m | 76 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 76 insertions(+)
 create mode 100644 sim/TurkeySimulation.m

diff --git a/sim/TurkeySimulation.m b/sim/TurkeySimulation.m
new file mode 100644
index 0000000..b696724
--- /dev/null
+++ b/sim/TurkeySimulation.m
@@ -0,0 +1,76 @@
+%for blog: take pictures of graphs of inside and outsite tempurature vs
+%time and oven temperature, also temperature vs radius of a properly cooked turkey.
+
+tic
+init=0;
+if init == 0
+clear all
+weight = 15/2.2;%lbs
+density = 800;
+% r0=.15;
+r0 = (weight/density/4*3/pi)^(1/3);
+deltar=r0/30;
+r=[deltar*1:deltar:r0];
+
+alpha = .412/2800/1000;
+h=9; %nusselt number
+k=.412+r*0;
+
+% weight = r0^3*4*pi/3*density;
+
+Tinf=600;
+Ti=20;
+T=Ti+r*0;
+waterlost=r*0;
+% T(end)=T(end)+1;
+end
+
+deltat = .1;
+time=.6*60*60/deltat;%minutes
+Tcenterlog=zeros(1,time);
+Tlog=Tcenterlog;
+
+for i = 1:time
+
+    dTtemp = (T(3:end)-T(1:end-2))./(2*deltar); %only works if r is linear
+    dTdr= [(T(2)-T(1))/deltar, dTtemp, (T(end)-T(end-1))/deltar];
+%     dTdr = diff(T)./diff(r);
+%     dTdr(end+1)=dTdr(end);
+    dTdr(end)= h*(Tinf-T(end))/k(end);
+    parenthesis=dTdr.*r.^2;
+    dPtemp = (parenthesis(3:end)-parenthesis(1:end-2))./(2*deltar); %again, linear only
+    dPdr = [(parenthesis(2)-parenthesis(1))/deltar, dPtemp, (parenthesis(end)-parenthesis(end-1))/deltar];
+    
+%     dPdr = diff(parenthesis)./diff(r);
+%     dPdr(end+1) = dPdr(end);
+    
+    
+% %     if kout > k(end)
+% %         k(end)=(kout-k(end))/3+k(end);
+% %     end
+%     kout=h*(Tinf-T(end))/dTdr(end);
+%     k(end)=kout;
+    
+    alpha = k/2800/density;
+    
+    dTdt=alpha.*dPdr./r.^2;
+    
+%     dTdt(1)=dTdt(1)/2;
+    
+    T=T+dTdt*deltat;
+    waterlost=waterlost+(10.^((T-20)/80)-1)/10;
+%     Tcenterlog(i)=T(1);
+%     Tlog(i)=T(end);
+
+end
+figure(1)
+plot(r,(T*1.8+32))
+
+figure(2)
+plot(r,waterlost)
+% cooktime=[1:time]*deltat/3600;
+% figure(2)
+% hold on
+% plot(cooktime,Tcenterlog*1.8+32)
+% plot(cooktime,Tlog*1.8+32)
+% toc
\ No newline at end of file