function = SunEarthMoon(years,framerate) %% SunEarthMoon % Final Project % By: Eric Keldrauk % Submission Date: 12/01/2008 % MATH 128A % % The purpose of this function is to plot the orbits of the sun, earth, and % moon over a specified time period (Note that time periods in excess of 20 % years will require large runtimes). % % Inputs: % years = Observation Time Period (years) % framerate = Video Speed (typically between 1 and 1000) %% Clean Up close all clc %% Initializaion x_earth = 147300000000; % [m] v_earth = 30257; % [m/s] r_sat = 384748000; % earth surface [m] r_earth = 6367000; % earth radius [m] v_sat = 1023; % relative velocity from earth [m/s] a = 5.145; % Angle to vertical (y) axis b = 90; % Angle to horizontal (x) axis in xz plane x_earth_o = [x_earth; 0; 0]; x_sun_o = [0; 0; 0]; x_sat_o = [x_earth+r_sat; 0; 0]; v_earth_o = [0; v_earth; 0]; v_sun_o = [0; 0; 0]; v_sat_o = v_sat*[cos(pi/180*b)*sin(pi/180*a); cos(pi/180*a); sin(pi/180*b)*sin(pi/180*a)] + v_earth_o; interval = years*[0 31536000]; %% Error Control h = [0.01 36000]; tol = 100000; Options.AbsTol = tol; Options.MaxStep = h(2); Options.InitialStep = h(1); %% Analysis ao = [x_earth_o; v_earth_o; x_sun_o; v_sun_o; x_sat_o; v_sat_o]; [t, x] = ode45(@earthfinal,interval,ao,Options); for i = 1:length(t) R1(i) = (x(i,13)-x(i,1)); R2(i) = (x(i,14)-x(i,2)); R3(i) = (x(i,15)-x(i,3)); R(i) = sqrt(R1(i)^2+R2(i)^2+R3(i)^2); end T_index_earth = find([1; x(:,4)].*[x(:,4);1]<=0); T_index_moon = find([1; R2(:)].*[R2(:); 1]<=0); for i = 4:length(T_index_earth) T_earth_semi(i-3) = (t(T_index_earth(i)-1)-t(T_index_earth(i-2)1))/24/60/60; end T_earth = mean(T_earth_semi); for i = 4:length(T_index_moon) T_moon_semi(i-3) = (t(T_index_moon(i)-1)-t(T_index_moon(i-2)-1))/24/60/60; end T_moon = mean(T_moon_semi); D_earth = 0; for i = 2:(T_index_earth(4)-1) D_earth = D_earth + sqrt((x(i,1)-x(i-1,1))^2+(x(i,2)-x(i-1,2))^2+(x(i,3)x(i-1,3))^2); end D_moon = 0; for i = 2:(T_index_moon(4)-1) D_moon = D_moon + sqrt((R1(i)-R1(i-1))^2+(R2(i)-R2(i-1))^2+(R3(i)-R3(i1))^2); end %% Plots q = framerate; scrsz = get(0,'ScreenSize'); figure('position', [0.05*scrsz(3) 0.05*scrsz(4) 0.75*scrsz(3) 0.85*scrsz(4)]) set(gcf,'name','Sun, Earth, and Moon Orbits') for i = 1:length(t)/q subplot(2,2,1) plot3(x(1:i*q,1),x(1:i*q,2),x(1:i*q,3),'g',x(1:i*q,7),x(1:i*q,8),x(1:i*q,9),'r ',x(1:i*q,13),x(1:i*q,14),x(1:i*q,15),'b') axis(1.1*[min(x(:,1)) max(x(:,1)) min(x(:,2)) max(x(:,2)) 2*min(x(:,15)) 2*max(x(:,15))]) xlabel('Universal X Coordinate (m)') ylabel('Universal Y Coordinate (m)') zlabel('Universal Z Coordinate (m)') title('Relative Orbits') legend('Earth','Sun','Moon') hold on plot3(x(i*q,1),x(i*q,2),x(i*q,3),'g-o',x(i*q,7),x(i*q,8),x(i*q,9),'ro',x(i*q,13),x(i*q,14),x(i*q,15),'b-o') hold off subplot(2,2,2) plot3(R1(1:i*q),R2(1:i*q),R3(1:i*q),'b',zeros(1,i*q),zeros(1,i*q),zeros(1,i*q) ,'g') axis(1.5*[min(R1) max(R1) min(R2) max(R2) min(R3) max(R3)]) xlabel('Universal X Coordinate (m)') ylabel('Universal Y Coordinate (m)') zlabel('Universal Z Coordinate (m)') title('Relative Moon Orbit About Earth') hold on plot3(R1(i*q),R2(i*q),R3(i*q),'b-o',0,0,0,'g-o') text(0,1.45*max(R2),1.40*max(R3),sprintf('Orbital Period, T = %3.5g days',T_moon)) text(0,1.45*max(R2),1.15*max(R3),sprintf('Orbital Circumference, D = %3.5g gigameters',D_moon*1e-9)) hold off subplot(2,2,3) plot(x(1:i*q,1),x(1:i*q,2),'g',x(1:i*q,7),x(1:i*q,8),'r') axis(1.5*[min(x(:,1)) max(x(:,1)) min(x(:,2)) max(x(:,2))]) xlabel('Universal X Coordinate (m)') ylabel('Universal Y Coordinate (m)') title('Relative Earth Orbit About Sun') hold on plot(x(i*q,1),x(i*q,2),'g-o',x(i*q,7),x(i*q,8),'r-o') text(1.45*min(x(:,1)),1.40*max(x(:,2)),sprintf('Orbital Period, T = %3.5g days',T_earth)) text(1.45*min(x(:,1)),1.25*max(x(:,2)),sprintf('Orbital Circumference, D = %3.5g gigameters',D_earth*1e-9)) text(1.45*min(x(:,1)),1.40*min(x(:,2)),sprintf('Time, t = %3.3g days',round(t(q*i)/24/60/60))) hold off subplot(2,2,4) plot(t(1:i*q)/(60*60*24),R(1:i*q)/1000,'b') axis([t(1)/24/60/60 t(end)/24/60/60 0.999*min(R)/1000 1.001*max(R)/1000]) xlabel('Time,t (days)') ylabel('Orbit Radius, R (km)') title('Moon-Earth Distance') hold on plot(t(i*q)/(60*60*24),R(i*q)/1000,'b-o') hold off end end drawnow %% Differential Equation Function function [udot]= earthfinal(t,u) m_earth = 5.9742e24; % [kg] m_sun = 1.98892e30; % [kg] m_sat = 11110; % [kg] G = 6.67300e-11; %[(m)^3(kg)^-1(s)^-2]; d_earth_sun = sqrt((u( 7,1)-u(1,1))^2+(u( 8,1)-u(2,1))^2+(u( 9,1)-u(3,1))^2); d_earth_sat = sqrt((u(13,1)-u(1,1))^2+(u(14,1)-u(2,1))^2+(u(15,1)-u(3,1))^2); d_sun_sat = sqrt((u(13,1)-u(7,1))^2+(u(14,1)-u(8,1))^2+(u(15,1)-u(9,1))^2); % Earth motion udot( 1,1) = u(4,1); udot( 2,1) = u(5,1); udot( 3,1) = u(6,1); udot( 4,1) = m_sun*G*(u(7,1)-u(1,1))/d_earth_sun^3 + m_sat*G*(u(13,1)u(1,1))/d_earth_sat^3; udot( 5,1) = m_sun*G*(u(8,1)-u(2,1))/d_earth_sun^3 + m_sat*G*(u(14,1)u(2,1))/d_earth_sat^3; udot( 6,1) = m_sun*G*(u(9,1)-u(3,1))/d_earth_sun^3 + m_sat*G*(u(15,1)u(3,1))/d_earth_sat^3; % Sun Motion udot( 7,1) = u(10,1); udot( 8,1) = u(11,1); udot( 9,1) = u(12,1); udot(10,1) = m_earth*G*(u(1,1)-u(7,1))/d_earth_sun^3 + m_sat*G*(u(13,1)u(7,1))/d_sun_sat^3; udot(11,1) = m_earth*G*(u(2,1)-u(8,1))/d_earth_sun^3 + m_sat*G*(u(14,1)u(8,1))/d_sun_sat^3; udot(12,1) = m_earth*G*(u(3,1)-u(9,1))/d_earth_sun^3 + m_sat*G*(u(15,1)u(9,1))/d_sun_sat^3; % Satellite Motion udot(13,1) = u(16,1); udot(14,1) = u(17,1); udot(15,1) = u(18,1); udot(16,1) = m_earth*G*(u(1,1)-u(13,1))/d_earth_sat^3 + m_sun*G*(u(7,1)u(13,1))/d_sun_sat^3; udot(17,1) = m_earth*G*(u(2,1)-u(14,1))/d_earth_sat^3 + m_sun*G*(u(8,1)u(14,1))/d_sun_sat^3; udot(18,1) = m_earth*G*(u(3,1)-u(15,1))/d_earth_sat^3 + m_sun*G*(u(9,1)u(15,1))/d_sun_sat^3; end %% HW4 Problem 1 TotalMinutes1 = 151; [Hours1,Minutes1,Message1] = myTimeConversion (TotalMinutes1) TotalMinutes2 = 84; [Hours2,Minutes2,Message2] = myTimeConversion (TotalMinutes2) TotalMinutes3 = 60; [Hours3,Minutes3,Message3] = myTimeConversion (TotalMinutes3) TotalMinutes4 = 181; [Hours4,Minutes4,Message4] = myTimeConversion (TotalMinutes4) TotalMinutes5 = 1; [Hours5,Minutes5,Message5] = myTimeConversion (TotalMinutes5) %% Problem 2 message = 'UC Berkeley is everything I hoped for.'; secretKey = 'abc'; ctext = Convert(message,secretKey,'encrypt') Convert(ctext,secretKey,'decrypt') %% Problem 3 temp1 = 50; unitsIn1 = 'F'; unitsOut11 = 'C'; unitsOut12 = 'K'; [newTemp11,msg11] = myTempConversion(temp1, unitsIn1, unitsOut11) [newTemp12,msg12] = myTempConversion(temp1, unitsIn1, unitsOut12) temp2 = 10; unitsIn2 = 'C'; unitsOut21 = 'F'; unitsOut22 = 'K'; [newTemp21,msg21] = myTempConversion(temp2, unitsIn2, unitsOut21) [newTemp22,msg22] = myTempConversion(temp2, unitsIn2, unitsOut22) temp3 = 30; unitsIn3 = 'K'; unitsOut31 = 'F'; unitsOut32 = 'C'; [newTemp31,msg31] = myTempConversion(temp3, unitsIn3, unitsOut31) [newTemp32,msg32] = myTempConversion(temp3, unitsIn3, unitsOut32) temp4 = 30; unitsIn4 = 'K'; unitsOut4 = 'K'; [newTemp4,msg4] = myTempConversion(temp4,unitsIn4,unitsOut4) %% Problem 4 a v = [1,0,3,4,5]; vpad1 = padrv(v,2) vpad2 = padrv(v,9) vpad3 = padrv(v,4) %% Problem 4 b b1 = unpad(vpad1) b2 = unpad(vpad2) b3 = unpad(vpad3) %% Problem 4 c p1 = [1,0,3]; p2 = [-1,2,1]; result9c1 = polyadd(p1,p2) p3 = [1,1,0,3]; p4 = [-1,2,1]; result9c2 = polyadd(p3,p4) %% Problem 5 f = @(x) cos(x.^2); df = @(x) -2*x.*sin(x.^2); h = 0.1; [df1est,df2est] = derivFCN(f,h); x = linspace(0,pi,100); plot(x,df(x),x,df1est(x),x,df2est(x)) legend('Exact Derivative','d1 Estimate','d2 Estimate'); %% Problem 6 f1h = @(z) cos(z)*sin(z); f2h = @(z) cos(z).*sin(z); L1 = 1; R1 = 2; N1 = 10; [sampleMean, sampleSdev, employV] = getStats(f1h,L1,R1,N1) [sampleMean, sampleSdev, employV] = getStats(f2h,L1,R1,N1) %% Problem 7 a A1 = [1.5,1.5,0]; A2 = [0,2.5,2.5]; A3 = [4,0,4]; B1 = myDiagnoser1(A1) B2 = myDiagnoser1(A2) B3 = myDiagnoser1(A3) %% Problem 7 b C1 D1 D2 D3 = = = = [1425,1432,1426]; C2 = [1432,1433,1419]; C3 = [1431,1432,1435]; myDiagnoser2(C1) myDiagnoser2(C2) myDiagnoser2(C3) %% HW4 Problem 1 TotalMinutes1 = 151; [Hours1,Minutes1,Message1] = myTimeConversion (TotalMinutes1) TotalMinutes2 = 84; [Hours2,Minutes2,Message2] = myTimeConversion (TotalMinutes2) TotalMinutes3 = 60; [Hours3,Minutes3,Message3] = myTimeConversion (TotalMinutes3) TotalMinutes4 = 181; [Hours4,Minutes4,Message4] = myTimeConversion (TotalMinutes4) TotalMinutes5 = 1; [Hours5,Minutes5,Message5] = myTimeConversion (TotalMinutes5) %% Problem 2 message = 'UC Berkeley is everything I hoped for.'; secretKey = 'abc'; ctext = Convert(message,secretKey,'encrypt') Convert(ctext,secretKey,'decrypt') %% Problem 3 temp1 = 50; unitsIn1 = 'F'; unitsOut11 = 'C'; unitsOut12 = 'K'; [newTemp11,msg11] = myTempConversion(temp1, unitsIn1, unitsOut11) [newTemp12,msg12] = myTempConversion(temp1, unitsIn1, unitsOut12) temp2 = 10; unitsIn2 = 'C'; unitsOut21 = 'F'; unitsOut22 = 'K'; [newTemp21,msg21] = myTempConversion(temp2, unitsIn2, unitsOut21) [newTemp22,msg22] = myTempConversion(temp2, unitsIn2, unitsOut22) temp3 = 30; unitsIn3 = 'K'; unitsOut31 = 'F'; unitsOut32 = 'C'; [newTemp31,msg31] = myTempConversion(temp3, unitsIn3, unitsOut31) [newTemp32,msg32] = myTempConversion(temp3, unitsIn3, unitsOut32) temp4 = 30; unitsIn4 = 'K'; unitsOut4 = 'K'; [newTemp4,msg4] = myTempConversion(temp4,unitsIn4,unitsOut4) %% Problem 4 a v = [1,0,3,4,5]; vpad1 = padrv(v,2) vpad2 = padrv(v,9) vpad3 = padrv(v,4) %% Problem 4 b b1 = unpad(vpad1) b2 = unpad(vpad2) b3 = unpad(vpad3) %% Problem 4 c p1 = [1,0,3]; p2 = [-1,2,1]; result9c1 = polyadd(p1,p2) p3 = [1,1,0,3]; p4 = [-1,2,1]; result9c2 = polyadd(p3,p4) %% Problem 5 f = @(x) cos(x.^2); df = @(x) -2*x.*sin(x.^2); h = 0.1; [df1est,df2est] = derivFCN(f,h); x = linspace(0,pi,100); plot(x,df(x),x,df1est(x),x,df2est(x)) legend('Exact Derivative','d1 Estimate','d2 Estimate'); %% Problem 6 f1h = @(z) cos(z)*sin(z); f2h = @(z) cos(z).*sin(z); L1 = 1; R1 = 2; N1 = 10; [sampleMean, sampleSdev, employV] = getStats(f1h,L1,R1,N1) [sampleMean, sampleSdev, employV] = getStats(f2h,L1,R1,N1) %% Problem 7 a A1 = [1.5,1.5,0]; A2 = [0,2.5,2.5]; A3 = [4,0,4]; B1 = myDiagnoser1(A1) B2 = myDiagnoser1(A2) B3 = myDiagnoser1(A3) %% Problem 7 b C1 D1 D2 D3 = = = = [1425,1432,1426]; C2 = [1432,1433,1419]; C3 = [1431,1432,1435]; myDiagnoser2(C1) myDiagnoser2(C2) myDiagnoser2(C3) @ This file can be opened as a Live Script. For more information, see greasing Live Scripts. 1 9696 Hi\" Problem 1 2 - TotalMinutesl = 151; 3 - [Hoursl,Minutesl,Message1] .=. myTimeConversion (TotalMinutesl) 4 - TotalMinutesZ = 84; 5 - [HoursZ,MinutesZ,Message2] .=. myTimeConversion (TotalMinutesZ) - TotalMinutesB = 60; - [Hours3,Minute53,Message3] .. myTimeConversion (TotaLMinutes3) - TotalMinutes4 = 181; [Hours4,Minutes4,Message4] .. myTimeConversion (TotaLMinutes4) ' :R 1- 0 Current Folder Command Window I Namel Error: Unbalanced or unexpected parenthesis or bracket. i_I $RECYCLE.BIN b i_I MathWorks >> Accept_now[1] E] 0 introduction.pdf Accept_now [1] m 1 simple random sampling.pdf r "f \"'\""" """'"" """""\" ' Error: Unbalanced or unexpected parenthesis or bracket. Details A >> Accept_now Undefined function or variable 'myTimeConversion'. Workspace E) Name 1 Value EEITotalMinutesl 151 Error in Accept now (line 3) g [Hoursl,Minutesl,Message1] = myTimeConversion (TotalMinutesl) %% HW4 Problem 1 TotalMinutes1 = 151; [Hours1,Minutes1,Message1] = myTimeConversion (TotalMinutes1) TotalMinutes2 = 84; [Hours2,Minutes2,Message2] = myTimeConversion (TotalMinutes2) TotalMinutes3 = 60; [Hours3,Minutes3,Message3] = myTimeConversion (TotalMinutes3) TotalMinutes4 = 181; [Hours4,Minutes4,Message4] = myTimeConversion (TotalMinutes4) TotalMinutes5 = 1; [Hours5,Minutes5,Message5] = myTimeConversion (TotalMinutes5) %% Problem 2 message = 'UC Berkeley is everything I hoped for.'; secretKey = 'abc'; ctext = Convert(message,secretKey,'encrypt') Convert(ctext,secretKey,'decrypt') %% Problem 3 temp1 = 50; unitsIn1 = 'F'; unitsOut11 = 'C'; unitsOut12 = 'K'; [newTemp11,msg11] = myTempConversion(temp1, unitsIn1, unitsOut11) [newTemp12,msg12] = myTempConversion(temp1, unitsIn1, unitsOut12) temp2 = 10; unitsIn2 = 'C'; unitsOut21 = 'F'; unitsOut22 = 'K'; [newTemp21,msg21] = myTempConversion(temp2, unitsIn2, unitsOut21) [newTemp22,msg22] = myTempConversion(temp2, unitsIn2, unitsOut22) temp3 = 30; unitsIn3 = 'K'; unitsOut31 = 'F'; unitsOut32 = 'C'; [newTemp31,msg31] = myTempConversion(temp3, unitsIn3, unitsOut31) [newTemp32,msg32] = myTempConversion(temp3, unitsIn3, unitsOut32) temp4 = 30; unitsIn4 = 'K'; unitsOut4 = 'K'; [newTemp4,msg4] = myTempConversion(temp4,unitsIn4,unitsOut4) %% Problem 4 a v = [1,0,3,4,5]; vpad1 = padrv(v,2) vpad2 = padrv(v,9) vpad3 = padrv(v,4) %% Problem 4 b b1 = unpad(vpad1) b2 = unpad(vpad2) b3 = unpad(vpad3) %% Problem 4 c p1 = [1,0,3]; p2 = [-1,2,1]; result9c1 = polyadd(p1,p2) p3 = [1,1,0,3]; p4 = [-1,2,1]; result9c2 = polyadd(p3,p4) %% Problem 5 f = @(x) cos(x.^2); df = @(x) -2*x.*sin(x.^2); h = 0.1; [df1est,df2est] = derivFCN(f,h); x = linspace(0,pi,100); plot(x,df(x),x,df1est(x),x,df2est(x)) legend('Exact Derivative','d1 Estimate','d2 Estimate'); %% Problem 6 f1h = @(z) cos(z)*sin(z); f2h = @(z) cos(z).*sin(z); L1 = 1; R1 = 2; N1 = 10; [sampleMean, sampleSdev, employV] = getStats(f1h,L1,R1,N1) [sampleMean, sampleSdev, employV] = getStats(f2h,L1,R1,N1) %% Problem 7 a A1 = [1.5,1.5,0]; A2 = [0,2.5,2.5]; A3 = [4,0,4]; B1 = myDiagnoser1(A1) B2 = myDiagnoser1(A2) B3 = myDiagnoser1(A3) %% Problem 7 b C1 D1 D2 D3 = = = = [1425,1432,1426]; C2 = [1432,1433,1419]; C3 = [1431,1432,1435]; myDiagnoser2(C1) myDiagnoser2(C2) myDiagnoser2(C3)