MATLAB太阳黑子活动周期.

太阳黑子的活动周期研究•班级：电控研14-2•姓名：庆春太阳黑子（太阳的暗黑斑点）•太阳黑子（sunspot）是在太阳的光球层上发生的一种太阳活动，是太阳活动中最基本、最明显的。一般认为，太阳黑子实际上是太阳表面一种炽热气体的巨大漩涡，温度大约为3000-4500℃。•因为其温度比太阳的光球层表面温度要低1000到2000摄氏度（光球层表面温度约为6000摄氏度），所以看上去像一些深暗色的斑点。太阳黑子很少单独活动，通常是成群出现。•黑子的活动周期为11年，活跃时会对地球的磁场产生影响，主要是使地球南北极和赤道的大气环流作经向流动，从而造成恶劣天气，使气候转冷。严重时会对各类电子产品和电器造成损害。•太阳黑子虽然颜色较深，但是在观测情况下，与太阳耀斑同样清晰显眼。天文学家把太阳黑子最多的年份称为“太阳活动峰年”，太阳黑子最少的年份称为“太阳活动谷年如何观测•一般的光学望远镜主要是观测太阳光球层的活动现象：黑子、耀斑和米粒组织。我们主要观测黑子和耀斑。•1．太阳黑子是光球上最明显的活动，在日面的形成、发展、消失、形态、多少、大小和分布等，都有一定的规律。黑子在日面上的东西分布是不对称的——东半球比西半球多，并且集中分布在日面中纬度区域。较完整的黑子是由本影和半影组成，中间深黑色的部分叫本影，四周暗淡的部分叫半影。另外大多数黑子都成群出现，多是沿东西方向分布着。一群中，西边较大的黑子叫前导黑子，东边较大的黑子叫后随黑子。就太阳整体而言，黑子活动有大约11年的周期变化。•2．光斑是在太阳光球层的东西边缘，往往出现一些比光球的背景更明亮的区域，这就是光斑。光斑常伴随日面边缘的黑子一起出现，也存在着11年的周期。•观测的方法：太阳观测有几种不同的方法。我们暂时只选用投影观测，以后进行照相观测。下面的内容中将重点介绍一下投影观测，简要说明照相观测。•投影观测•一、测前准备•（1）观测纸：印有日本威信公司生产的太阳观测网格图，观测者根据当日的B0、P值选择，共有八种（在下文中介绍）。•（2）描图纸：印有与观测纸网格图同样大小稀疏网格图（直径为10厘米，B0=0度，P=0度），也可以用白而薄的纸（纸上有直径10厘米的圆）•（3）天顶棱镜：若使用的太阳投影板平面与镜筒平行，需使用天顶棱镜改变光路。•（4）目镜：不能使用胶合目镜，用惠更斯目镜较为合适。不能使用高倍目镜。•（5）滤光片：ND中性滤光片，有ND2、ND4、ND8多种。•（6）天文普及年历：观测当天的B0、L0、P值都需从当年的天文年历中查得。•（7）常用的工具：铅笔、手表、红黄彩笔、尺、圆规、投影仪、小磁铁或夹子。•（8）带上前一天的描点图。•1、将望远镜极轴大致对准北极，镜筒指向太阳，观察镜筒的影子，当影子最小时，说明太阳已经大致对好。•！！！请注意，千万不能用眼睛通过望远镜直接寻找太阳，寻星镜也最好不要加装。•2、记录透明度、宁静度、云量、描图开始时间等观测背景数据。关于宁静度和透明度的规定是这样的：范围是0-5最好是5，最差是0，一般都要在2.5以上才可以观测。宁静度与风有关，一般4级以下可以定为3。透明度与污染等有关，取值看经验了。•3、装上、目镜、投影板，调好焦距，将描图纸固定在投影板上（我们用小的吸铁石或夹子），上下移动板，使太阳的像与描图纸上的圆面等大。•4、不要跟踪太阳，选择一颗小而清楚的黑子，观察其在描图纸上的视运动轨迹，确定东西线。•5、若使用网格图来描点，只需旋转描图纸，使黑子的视运动轨迹方向（太阳的东西线方向）与描图纸中的纬线平行。使用白纸只需在数据分析时画出东西线即可。•6、确定日面坐标方向•①若投影板与镜筒平行，由于描图纸上的像是经过折射镜和棱镜两次改变光路得到的，故太阳黑子视运动方向恰与日面坐标中的东西走向大体一致，再按“上北下南，左西右东”，日面坐标方向就确定了。•②若投影板与镜筒垂直，由于只经过一次改变光路，方向完全相反：按上南下北，左东右西。即可确定日面坐标方向。•7、把投影相和描图纸的圆对齐并开始跟踪，在日面上找一颗形状较为规则的大黑子，适当调整望远镜的赤经、赤纬旋钮，当投影日面与描图纸上的圆重合时迅速描下该黑子的位置，其他黑子就可以以它为主参照黑子定位。另外边看边记黑子群数、个数、形态、位置等。此时要对照前一天的黑子描点图，并确定出黑子群的类型（采用苏黎式分类法）•8、描图时要迅速，手不要压屏，头不要碰目镜。•9、其他黑子的描图需按以下步骤进行：•①黑子需自西向东描。•②先描黑子半影，后描本影。•③先描前导黑子，后描后随黑子。•10、描出光斑。•11、记录描图结束的时间，以描图开始结束时间的中间时刻，记为描图时间。数据•见附表•太阳黑子的活动呈现周期性变化的，沃尔夫(R.Wolfer)根据在过去的288年(1700年~1987年)间每年太阳黑子出现的数量和大小的观测数据推算出11年的周期规律。数据处理•离散时间信号的傅里叶变换•DiscreteTimeFourierTransform,DTFT•DFS是周期序列的离散傅里叶级数•DTFT是非周期序列的傅里叶变换，称离散时间傅里叶变换，其频谱是连续的函数•DFT是有限长序列的离散傅里叶变换，是对其DTFT的等间隔抽样，是离散的频谱•DFT是DFS的主值序列，是非周期的。而DFS是DTFT的频域内的抽样。•FFT是DFT的一种高效快速算法，也称作快速傅里叶变换。离散傅立叶变换（DFT）•公式•这一对式子，左、右两边都是离散的，有限长，因此可方便地用来实现频谱分析。•但使用时，一定要想到，它们均来自DFS,即•x（n）和X（k）都是周期的！•FFT是一种DFT的高效算法，称为快速傅立叶变换（fastFouriertransform）。FFT算法可分为按时间抽取算法和按频率抽取算法，先简要介绍FFT的基本原理。从DFT运算开始，说明FFT的基本原理。•在该实验中，对Wolfer数序列做FFT变换后得到Y（长度为n），只取其前n/2个数据的功率谱密度的估计值|Y|^2。原因是时域为离散的实序列的傅立叶变换对应于具有周期性且偶对称的频域特性，因此Y的前n/2个数据已经包含了Wolfer数的全部信息。根据DFT的频域单位ω=2π/N*k与DTFT的频域单位ω=2𝜋𝑓𝑠𝑓的表达关系式以及k与f对应关系，可以看出k与f呈线性关系；同样地，因为Y的前n/2个数据已经包含了Wolfer数的全部信息，只取前n/2个数据分析功率-频率图时，对应的横坐标时应取。（注：𝑓𝑠为采样频率）实现步骤及内容•(1)读取数据并绘制Wolfer图，观察太阳黑子的活动现象。程序如下：•loadsunsb.dat%读取数据•year=sunb(:,1);%读取年份信息•wolfer=sunb(:,2);%读取黑子活动数据•figure%新建图像•plot(year,wolfer)%画出时域图•xlabel('Years');ylabel('SunbData');title('SunbData')%标注横纵坐标和标题•pause•为清楚起见，画出部分时间段(如实验中采取最近50年)的Wolfer图进行观察。程序如下：•figure•plot(year(261:311),wolfer(261:311),'b.-');%用进50年数据做图•xlabel('Years');ylabel('SunbData');title('Atthelast50years')%标注横纵坐标和标题(2)对已经得到的Wolfer数应用FFT技术分析它的变化规律。程序如下：Y=fft(wolfer);%对全部数据做FFTY(1)=[];%由于算法问题，Y(1)为所有fft数值之和，因此需舍弃第一个点•pause观察Wolfer数在FFT分析后得到的在复平面中的分布图。程序如下：figureplot(Y,'ro')%在复平面做图，空心点title('FourierCoefficientsintheComplexPlane');xlabel('RealAxis');%标注横坐标为’实部’ylabel('ImaginaryAxis');%标注纵坐标为’虚部’pause(3)功率与频率的关系曲线则被定义为周期图。用周期图法对Wolfer数进行功率谱密度分析。程序：figuren=length(Y);%取FFT结果长度为npower=abs(Y(1:n/2)).^2;%取前n/2个数据求其功率nyquist=1/2;%取最大频率为0.5freq=(1:n/2)/(n/2)*nyquist;%将FFT转换为DTFT%假设采样频率为Fs，采样点数为N，做FFT之%后，某一点n（n从1开始）表示的频率为：%Fn=(n-1)*Fs/N；该点的模值除以N/2就是对应。真实频率plot(freq,power)%画周期图xlabel('cycles/year');title('Periodogram')pause为清楚起见，对功率和频率的前50个分量作它的周期图。程序如下：figureplot(freq(1:50),power(1:50))%为观察方便xlabel('cycles/year')pause(4)确定出太阳黑子的活动周期。为清楚起见，画出功率与周期（频率的倒数）的关系曲线图。程序如下：figureperiod=1./freq;%将横坐标单位转换为年/周期plot(period,power);%画周期图axis([05002e+7]);ylabel('Power');xlabel('Period(Years/Cycle)');%单位为年/周期pause(5)在功率与周期关系曲线图中标出功率的最高点，该位置对应的周期则为太阳黑子活动周期。holdon;index=find(power==max(power));%找到频率最大点，该点横坐标即为太阳黑子周期mainPeriodStr=num2str(period(index));plot(period(index),power(index),'r.','MarkerSize',25);%用实心点指出该点text(period(index)+2,power(index),['Period=',mainPeriodStr]);%文字标注该点holdoff;谢谢