附录 1：附件 1-3 实际影长
附件 1 附件 2 附件 3
北京时间 影长 北京时间 影长 北京时间 影长
14:42 1.14962583 12:41 1.247256 13:09 3.533142 14:45 1.18219898 12:44 1.222795 13:12 3.546768 14:48 1.21529696 12:47 1.198921 13:15 3.561798 14:51 1.24905105 12:50 1.175429 13:18 3.578101 14:54 1.28319534 12:53 1.15244 13:21 3.595751 14:57 1.31799315 12:56 1.129917 13:24 3.614934 15:00 1.35336405 12:59 1.107835 13:27 3.635426 15:03 1.38938709 13:02 1.086254 13:30 3.657218 15:06 1.42615286 13:05 1.065081 13:33 3.680541 15:09 1.46339985 13:08 1.044446 13:36 3.705168 15:12 1.50148162 13:11 1.024264 13:39 3.731278 15:15 1.54023182 13:14 1.00464 13:42 3.758918 15:18 1.57985332 13:17 0.985491 13:45 3.788088 15:21 1.62014452 13:20 0.96679 13:48 3.818701 15:24 1.66127061 13:23 0.948585 13:51 3.85081 15:27 1.70329063 13:26 0.930928 13:54 3.884585 15:30 1.74620591 13:29 0.913752 13:57 3.919912 15:33 1.79005092 13:32 0.897109 14:00 3.956876 15:36 1.83501427 13:35 0.880974 14:03 3.995535 15:39 1.880875 13:38 0.865492 14:06 4.035751 15:42 1.92791845 13:41 0.850504 14:09 4.077863 附录 2：角度变化
角度变化
0.455535 0.343812
0.440947 0.330544
0.424696 0.326406
0.413645 0.316915
0.398599 0.311967
0.391943 0.306906
0.37772 0.298737
0.365601 0.292755
0.358179 0.287602
0.34811 0.28526
附录三 clear;
clc;
t = sym('t');
fai=[39,54,26];
seita=[116,23,29];
date='10/22';
h=3;
%% 日期的处理及赤纬角的计算 d=datenum(date,'mm/dd');
d0=datenum('1/1','mm/dd');
n=d-d0+1; % n 为一年中日期的序号 delta=23.45*sind(360*(284+n)/365);
delta=delta*pi/180;
%% 时间的处理及时角的计算 seita=dms2degrees(seita);
t0=(t-(120-seita)*4/60); % t0 为当地时间 omega=15*(abs(t0)-12); % omega 为时角 omega=omega*pi/180;
%% 求解高度角 fai=dms2degrees(fai);
fai=fai*pi/180;
hd=sin(fai)*sin(delta)+cos(fai)*cos(delta)*cos(omega);
alfa=asin(hd);
L=h/tan(alfa);
%% 求出最低点的坐标 fun=matlabFunction(L,'vars',t);
[X,Y]=fminsearch(fun,9,15);
fprintf('最短的影长为%d\n',Y);
%% 画出太阳影子变化曲线
%subplot(1,2,1);
ezplot(L,[9,15]);
title('天安门广场 9 点到 15 点太阳影长变化曲线');
xlabel ('北京时间');
ylabel ('影子的长度');
hold on
axis([9 15 3 8])
% %% 误差检验
% x=[7,7.5,8,8.5,9,9.5,10,10.5,11,11.5,12,12.5,13,13.5,14,14.5,15,15.5,16,16.5,17];
%
y=[38.75,17.35,11.21,8.29,6.59,5.5,4.76,4.25,3.92,3.73,3.67,3.74,3.94,4.28,4.8,5.56,6 .68,8.43,11.47,17.99,42.12];
% plot(x,y-0.4,'r-'); %画出数据图
% plot(x,y+1.2,'r-'); %画出数据图
% legend('MATLAB 计算得到的曲线','误差域');
% %% 画出太阳影子变化率曲线
% subplot(1,2,2);
% l=diff(L,t);
% ezplot(l,[9,15]);
% title('天安门广场太阳影长 9 点到 15 点变化率曲线');
% xlabel ('北京时间');
% ylabel ('影长变化率');
% xi=9:0.4:15;
% t=9:15;
% yi=interp1(t,eval(L),xi, 'spline');
% alfa=atan(h./yi);
% alfa=alfa.*180./pi;
附录四 clear;
clc;
%% 求出函数表达式
t = sym('t'); % t 为北京时间 fai = sym('fai'); % fai 为纬度 seita = sym('seita'); % seita 为经度
date = sym('date'); % date 为一年中日期的序号 h= sym('h'); % h 为杆长
% 时间的处理及时角的计算
t0=(t-(120-seita)*4/60); % t0 为当地时间 omega=15*(abs(t0)-12); % omega 为时角 omega=omega*pi/180;
% 日期的处理及赤纬角的计算 delta=23.45*sin(2*pi*(284+date)/365);
delta=delta*pi/180;
% 求解高度角 fai=fai*pi/180;
hd=sin(fai)*sin(delta)+cos(fai)*cos(delta)*cos(omega);
alfa=asin(hd);
global L;
L=h/tan(alfa)
%% 求出太阳影长与日期的的关系 subplot(2,1,1);
t = 12; % t 为北京时间 fai = 39; % fai 为纬度 seita = 116; % seita 为经度
date = sym('date'); % date 为一年中日期的序号 h= 3; % h 为杆长
L1=subs(L);
ezplot(L1,[1,365]);
title('太阳影长与日期的关系');
xlabel ('日期序号');
ylabel ('影子的长度');
l1=diff(L1, date);
ezplot(l1,[1,365]);
title('日期对太阳影长的影响率');
xlabel ('日期的变化');
ylabel ('太阳影长的变化率');
date=1:365;
lv=eval(l1);
d=['03/20';'03/21';'06/21';'06/22';'09/22';'09/23';'12/21';'12/22'];
d0=datenum('1/1','mm/dd');
for i=1:size(d,1)
day(i)=datenum(d(i,:),'mm/dd');
n=day-d0+1; % n 为一年中日期的序号 end
dlv=lv(n);
%% 求出太阳影长与直杆长度的关系 figure(2);
subplot(2,2,1);
t = 12; % t 为北京时间 fai = 39; % fai 为纬度 seita = 116; % seita 为经度
date = 295; % date 为一年中日期的序号 h= sym('h'); % h 为杆长
L2=subs(L);
ezplot(L2,[0,10]);
title('太阳影长与直杆长度的关系');
xlabel ('直杆的长度');
ylabel ('影子的长度');
% 直杆长度的变化对太阳影长的影响率 figure(3);
subplot(2,2,1);
l2=diff(L2, h);
ezplot(l2);
title('直杆长度的变化对太阳影长的影响率');
xlabel ('直杆的长度的变化长度');
ylabel ('太阳影长的变化率');
%% 求出太阳影长与经度的的关系 figure(2);
subplot(2,2,2);
t = 12; % t 为北京时间 fai = 39; % fai 为纬度 seita = sym('seita'); % seita 为经度
date = 295; % date 为一年中日期的序号 h= 3; % h 为杆长
L3=subs(L);
ezplot(L3,[60,180]);
title('太阳影长与经度的关系');
xlabel ('经度:60°E- 180°E');
ylabel ('影子的长度');
% 经度对太阳影长的影响率 figure(3);
subplot(2,2,2);
l3=diff(L3, seita);
ezplot(l3,[60,180]);
title('经度对太阳影长的影响率');
xlabel ('经度的变化');
ylabel ('太阳影长的变化率');
%% 求出太阳影长与纬度的的关系 figure(2);
subplot(2,2,3);
t = 12; % t 为北京时间 fai = sym('fai'); % fai 为纬度 seita = 116; % seita 为经度
date = 295; % date 为一年中日期的序号 h= 3; % h 为杆长
L4=subs(L);
ezplot(L4,[-30,30]);
title('太阳影长与纬度的关系');
xlabel ('纬度：30°S- 30°N');
ylabel ('影子的长度');
% 纬度对太阳影长的影响率 figure(3);
subplot(2,2,3);
l4=diff(L4, fai);
ezplot(l4,[-30,30]);
title('纬度对太阳影长的影响率');
xlabel ('纬度的变化');
ylabel ('太阳影长的变化率');
%% 求出太阳影长与时间的的关系 figure(2);
subplot(2,2,4);
t = sym('t'); % t 为北京时间 fai = 39; % fai 为纬度 seita = 116; % seita 为经度
date = 295; % date 为一年中日期的序号 h= 3; % h 为杆长
L5=subs(L);
ezplot(L5,[8.5,16.5]);
title('太阳影长与时间的关系');
xlabel ('一天内的时刻');
ylabel ('影子的长度');
% 一天中时间变化对太阳影长的影响率 figure(3);
subplot(2,2,4);
l5=diff(L5, t);
ezplot(l5,[8.5,16.5]);
title('一天中时间变化对太阳影长的影响率');
xlabel ('一天中时间变化');
ylabel ('太阳影长的变化率');
附录五 clear;
clc;
[a,b,c]=xlsread('附件 1-3.xls','附件 1','A4:C24');
[~,~,n]=xlsread('附件 1-3.xls','附件 1','A2');
%% 根据数据拟合出曲线 x=a(:,1);
y=a(:,2);
plot(x,y);
p = polyfit(x,y,2);
y1= polyval(p,x);
plot(x,y,'*r',x,y1,'-b');
legend('题目所给影端的位置','拟合得到的影端曲线');
title('太阳影长变化曲线');
xlabel ('横坐标 X(米)');
ylabel ('纵坐标 Y(米)');
%% 计算出影长 l=sqrt(x.^2+y.^2);
for i=1:length(c)
T{i}=datevec(c{i,1},'HH:MM:SS');
t(i,1)=T{i}(4)+T{i}(5)/60+T{i}(6)/3600;
end
%% 日期的处理
d=datenum(n,'yy/mm/dd');
d0=datenum('2015/1/1','yy/mm/dd');
date=d-d0+1; % n 为一年中日期的序号
%% 寻找经纬度 m=1;
CAL=zeros(21,2);
e0=zeros(1,2);
H=zeros(1,2);
FAI=zeros(1,2);
SEITA=zeros(1,2);
for h=1:10
for fai=0:0.6:60
for seita=15:0.6:150
cal = (h.*(1 - (sin((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600)...
.*sin((pi.*fai)./180) + cos((pi.*fai)./180).*...
cos((pi.*(15.*abs(seita./15 + t - 8) - 180))./180).*...
cos((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600)).^2).^(1./2))./...
(sin((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600).*sin((pi.*fai)./180) +...
cos((pi.*fai)./180).*cos((pi.*(15.*abs(seita./15 + t - 8) - 180))./180).*...
cos((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600));
error = sum(abs(cal-l)./l)/length(l);
if error<0.002 CAL(:,m)=cal;
e0(m)=error;
H(m)=h;
FAI(m)=fai;
SEITA(m)=seita;
m=m+1;
end end end end
%% 检验所求得经纬度下影长与实际影长的关系 global L;
% 第一个点 T=t;
figure(2);
subplot(1,2,1) seita=SEITA(1);
fai=FAI(1);
h=H(1);
t=sym('t');
L1=subs(L);
for i=1:21 t=T(i);
l1(i)=eval(L1);
end
stem(T,l1-l')
title('A 地所求影子长度的差值');
xlabel ('北京时间');
ylabel ('差值');
% ezplot(L1,[T(1),T(21)]);
% title('东经 109.2°北纬 19.2°点太阳影长与北京时间的关系');
% xlabel ('北京时间');
% ylabel ('影子的长度');
% hold on
% plot(T,l,'r');
% legend('所求得地点的影长与时间的关系曲线','题目数据拟合出的影长与时间 曲线',0);
% t=T;
% xi=14.7:0.05:15.7;
% yi=interp1(t,eval(L),xi, 'spline');
% alfa=atan(h./yi);
% alfa=alfa.*180./pi;
% 第二个点
seita=SEITA(2);
% yi=interp1(t,eval(L),xi, 'spline');
% alfa=atan(h./yi);
'8:10';
'8:00'];
n='1990/5/8';
T=datevec(t,'HH:MM');
t=T(:,4)+T(:,5)./60+T(:,6)./3600;
d=datenum(n,'yy/mm/dd');
d0=datenum('1990/1/1','yy/mm/dd');
date=d-d0+1; % n 为一年中日期的序号
%% 寻找经纬度 e0=inf;
for h=1:0.4:20 for fai=0:60
for seita=15:150
cal = (h.*(1 - (sin((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600)...
.*sin((pi.*fai)./180) + cos((pi.*fai)./180).*...
cos((pi.*(15.*abs(seita./15 + t - 8) - 180))./180).*...
cos((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600)).^2).^(1./2))./...
(sin((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600).*sin((pi.*fai)./180) +...
cos((pi.*fai)./180).*cos((pi.*(15.*abs(seita./15 + t - 8) - 180))./180).*...
cos((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600));
error = sum(abs(cal-l)./l)/length(l);
if error<e0 CAL=cal;
e0=error;
H=h;
FAI=fai;
SEITA=seita;
DATE=date;
end end end end 附录八 附录九 clear;
clc;
[a,b,c]=xlsread('附件 1-3.xls','附件 2','A4:C24');
%% 根据数据拟合出曲线
% 附件二 x1=a(:,1);
y1=a(:,2);
plot(x1,y1);
p = polyfit(x1,y1,2);
y1= polyval(p,x1);
plot(x1,y1,'*r',x1,y1,'-b');
title('太阳影长变化曲线');
xlabel ('横坐标 X(米)');
ylabel ('纵坐标 Y(米)');
% 计算出影长与时间 l1=sqrt(x1.^2+y1.^2);
for i=1:length(c)
T1{i}=datevec(c{i,1},'HH:MM:SS');
t1(i,1)=T1{i}(4)+T1{i}(5)/60+T1{i}(6)/3600;
end T1=t1;
figure(2);
% 时间与影长的拟合 p=polyfit(t1,l1,2);
y1=polyval(p,t1);%命令 polyval 多项式函数的预测值 t1=sym('t1');
z1=eval(sym('0.0982*t1^2-2.9856*t1+23.3233'));
% 求出最低点的坐标
fun=matlabFunction(z1,'vars',t1);
[X1,Y1]=fminsearch(fun,7,22);
s=120-15*(X1-12);
t=T1;
e0=inf;
for h=1:5
for fai=0:60
for date=1:365 for seita=s:s+10
cal = (h.*(1 - (sin((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600)...
.*sin((pi.*fai)./180) + cos((pi.*fai)./180).*...
cos((pi.*(15.*abs(seita./15 + t - 8) - 180))./180).*...
cos((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600)).^2).^(1./2))./...
(sin((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600).*sin((pi.*fai)./180) +...
cos((pi.*fai)./180).*cos((pi.*(15.*abs(seita./15 + t - 8) - 180))./180).*...
cos((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600));
error = sum(abs(cal-l1)./l1)/length(l1);
if error<e0 CAL=cal;
e0=error;
H=h;
FAI=fai;
DATE=date;
SEITA=seita;
end end end end
end global L;
seita=SEITA;
date=DATE;
fai=FAI;
h=H;
t=sym('t');
L1=subs(L);
ezplot(L1,[T1(1),T1(21)]);
title('东经 78.9°北纬 41°点太阳影长与北京时间的关系');
xlabel ('北京时间');
ylabel ('影子的长度');
hold on plot(T1,l1,'r');
clear;
clc;
[d,e,f]=xlsread('附件 1-3.xls','附件 3','A4:C24');
%% 根据数据拟合出曲线
% 附件二 x2=d(:,1);
y2=d(:,2);
plot(x2,y2);
p = polyfit(x2,y2,2);
y2= polyval(p,x2);
plot(x2,y2,'*r',x2,y2,'-b');
legend('题目所给影端的位置','拟合得到的影端曲线');
title('太阳影长变化曲线');
xlabel ('横坐标 X(米)');
ylabel ('纵坐标 Y(米)');
% 计算出影长与时间 l2=sqrt(x2.^2+y2.^2);
for i=1:length(f)
T2{i}=datevec(f{i,1},'HH:MM:SS');
t2(i,1)=T2{i}(4)+T2{i}(5)/60+T2{i}(6)/3600;
end T2=t2;
figure(2);
% 时间与影长的拟合 p=polyfit(t2,l2,2);
y2=polyval(p,t2);%命令 polyval 多项式函数的预测值 t2=sym('t2');
z2=eval(sym('0.2962*t2^2-7.5445*t2+51.5213'));
% 求出最低点的坐标
fun=matlabFunction(z2,'vars',t2);
[X2,Y2]=fminsearch(fun,7,22);
s=120-15*(X2-12);
t=T2;
for h=1:5
for fai=0:60
for date=1:365 for seita=s:s+10
cal = (h.*(1 - (sin((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600)...
.*sin((pi.*fai)./180) + cos((pi.*fai)./180).*...
cos((pi.*(15.*abs(seita./15 + t - 8) - 180))./180).*...
cos((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600)).^2).^(1./2))./...
(sin((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600).*sin((pi.*fai)./180) +...
cos((pi.*fai)./180).*cos((pi.*(15.*abs(seita./15 + t - 8) - 180))./180).*...
cos((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600));
error = sum(abs(cal-l2)./l2)/length(l2);
if error<e0 CAL=cal;
e0=error;
H=h;
FAI=fai;
DATE=date;
SEITA=seita;
end end end end end global L;
seita=SEITA;
date=DATE;
fai=FAI;
h=H;
t=sym('t');
L2=subs(L);
ezplot(L2,[T2(1),T2(21)]);
title('东经 109.9°北纬 28°点太阳影长与北京时间的关系');
xlabel ('北京时间');
ylabel ('影子的长度');
hold on plot(T2,l2,'r');
附录十 clear;
clc;
fileName = 'Appendix4.avi';
obj = VideoReader(fileName);
%% 从视频中提取出图片并计算出端点的像素坐标 n=0;
for i=1:3000:obj.NumberOfFrames
n=n+1;
video{n} = read(obj, i); % first frame only 获取视频中的 21 帧图片 %figure,I=imshow(video{n});
%saveas(I,strcat(num2str(n),'.jpg'));
end for i=1:n
x0(i,:) = point(video{i});
end
%% 在图片上描出几个时刻影子的端点
last = read(obj, inf); % last frame only 获取视频中的最后一帧图片 b1 = last(810:930,800:1700,:);
subplot(5,1,1),imshow(b1);
title('彩色图片');
c1 = rgb2gray(b1);
subplot(5,1,2),imshow(c1);
title('灰度图片');
d1 = im2bw(c1,0.81);
subplot(5,1,3),imshow(d1);
title('二值图片');
[z,l] = bwlabel(~d1);
index = z==z(50,20);
d1(~index) = 1;
subplot(5,1,4),imshow(d1);
title('二值降噪图片');
subplot(5,1,5),imshow(d1);
title('影端点的轨迹');
hold on,plot(x0(:,2),x0(:,1),'r*');
%% 选取的点的时间 t=[ '8:54:06';
'8:56:06';
'8:58:06';
'9:00:06';
'9:02:06';
'9:04:06';
'9:06:06';
'9:08:06';
'9:10:06';
'9:12:06';
'9:14:06';
'9:16:06';
'9:18:06';
'9:20:06';
'9:22:06';
'9:24:06';
'9:26:06';
'9:28:06';
'9:30:06';
'9:32:06';
];
T=datevec(t,'HH:MM:SS');
t=T(:,4)+T(:,5)./60+T(:,6)./3600;
%% 对图中直杆的长度处理 h=2;
b2 = last(200:900,760:980,:);
%figure,imshow(b2) c2 = rgb2gray(b2);
figure,imshow(c2) hold on
x1=[133,8];
x2=[110,673];
x3=[127,273];
x4=[114,273];
plot(x1(1),x1(2),'r*');
plot(x2(1),x2(2),'y*');
plot(x3(1),x3(2),'b*');
plot(x4(1),x4(2),'g*');
g = sqrt(sum((x1-x3).^2)) + sqrt(sum((x2- x4).^2));
y1=[56,686];
y2=[158,686];
dz1 = sqrt(sum((y1-y2).^2));
dz=h/g*dz1;
%% 利用杆座进行坐标变换 B=[ 14,77;
121,76;
118,56;
19,57;
];
A=[ 14,75;
121,75;
121,-32;
14,-32;
];
TForm = cp2tform(B,A,'projective');
X=tformfwd(TForm,[x0(:,2),x0(:,1)]);
dz2 = sqrt(sum((A(1,:)-A(2,:)).^2));
zb=X.*(dz/dz2);
%% 根据数据拟合出曲线 figure(3);
x=zb(:,1);
y=zb(:,2);
plot(x,y);
p = polyfit(x,y,2);
y1= polyval(p,x);
plot(x,y,'*r',x,y1,'-b');
legend('图片提取得到的影端位置','拟合得到的影端曲线');
title('太阳影长变化曲线');
xlabel ('横坐标 X(米)');
ylabel ('纵坐标 Y(米)');
%% 计算出影长 l=sqrt(x.^2+y.^2);
%% 日期的处理 n='2015/7/13';
d=datenum(n,'yy/mm/dd');
d0=datenum('2015/1/1','yy/mm/dd');
date=d-d0+1; % n 为一年中日期的序号
%% 寻找经纬度
% m=1;
% CAL=zeros(21,2);
% e0=zeros(1,2);
% H=zeros(1,2);
% FAI=zeros(1,2);
% SEITA=zeros(1,2);
e0=inf;
for fai=0:60
for seita=0:180
cal = (h.*(1 - (sin((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600)...
.*sin((pi.*fai)./180) + cos((pi.*fai)./180).*...
cos((pi.*(15.*abs(seita./15 + t - 8) - 180))./180).*...
cos((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600)).^2).^(1./2))./...
(sin((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600).*sin((pi.*fai)./180) +...
cos((pi.*fai)./180).*cos((pi.*(15.*abs(seita./15 + t - 8) - 180))./180).*...
cos((469.*pi.*sin((2.*pi.*(date + 284))./365))./3600));
error = sum(abs(cal-l)./l)/length(l);
if error<e0 e0=error;
CAL=cal;
H=h;
FAI=fai;
SEITA=seita;
end end end global L;
T=t;
figure(4);
seita=SEITA;
fai=FAI;
h=H;
t=sym('t');
L1=subs(L);
ezplot(L1,[T(1),T(21)]);
hold on plot(T,l,'r');
title('太阳影长与北京时间的关系');
ylabel ('影子的长度');
legend('所求得地点的影长与时间的关系曲线','图片提取得到的影长与时间曲线 ',0);
gg=[1.6;1.7;1.7;1.7;1.8;1.8;1.8;1.9;1.9;2;2.2;2.2;2.3];
wd=[43;39;40;41;38;39;40;39;40;40;36;37;39];
jd=[101;103;103;103;105;105;105;107;107;109;112;112;114];
cz=sqrt((wd-fai).^2+(jd-seita).^2);
figure(5);
gc=gg-h;
%plot(gc,cz)
p = polyfit(gc,cz,2);
y1= polyval(p,gc);
plot(gc,y1,'-b');
title('直杆高度变化对经纬度的影响');
xlabel ('直杆的高度差');
ylabel ('经纬度的误差');