2.Delimiter Matching

Prerequisites, Goals, and Outcomes

Prerequisites: Students should have mastered the following prerequisite skills. 

· Stack - Understanding of the implementation of stack. 

· Templates - Creation and use of template classes

· File IO - Basic file IO skills 

Goals: This assignment is designed to reinforce the student's understanding of the structure and use of stack as containers. 

Outcomes: Students successfully completing this assignment would master the following outcomes. 

· Familiarize how to implement a template class 

· Familiarize the structure and the use of the stack 

Background

Compilers check your programs for syntax errors, but frequently a lack of one symbol (such as a missing brace or comment starter) will cause the compiler to spill out a hundred lines of diagnostics without identifying the real error. 

A useful tool in this situation is a program that checks whether everything is balanced. Thus, every right brace, bracket, and parenthesis must correspond to their left counterparts. The sequence [()] is legal, but [(]) is wrong. Obviously, it is not worthwhile writing a huge program for this, but it turns out that it is easy to check these things. 

Description

Make an empty stack. Read characters until end of file. If the character is an open anything, push it onto the stack. If it is a close anything, then if the stack is empty report an error. Otherwise, pop the stack. If the symbol popped is not the corresponding opening symbol, then report an error. At end of file, if the stack is not empty report an error. 

Besides the balancing of the parentheses, brackets, and bracesthe Annotation (/**/) should be balanced. When a /* is found, all the letters behind will be ignored until a */ is found. If it is not found at the end of the file, an error will be reported.

Some of the contents of the stack are shown in the figure below.

Tasks

To complete this assessment, you need to define the implementation of class Stack and complete the main.cpp file in which the delimiter matching is defined. 

Following is an ordered list of steps that serves as a guide to completing this assessment. Work and test incrementally. Save often.

1. Begin by the creation of the class StackYou can implement a stack through either an array or a list. But this class should be defined as a template class. And it should contain the following basic operation. 

void Pop():  Pop an element from the top of the stack 

void Push(const T & e):  Push this element into the stack.

T & Top(): Return the top element of the stack.

bool Isempty() const: If the stack is empty, return true; otherwise return false.

int Size() const:  Return the number of the elements in the stack.

2. Nextcreate the file main.cpp and complete the implementation of delimiter matching. The detail steps are as follows.

Get a file name from the command line parameter, which specifies the file name needing delimiter matching. Invoke the delimiterMatching function using the file name.

Complete the function delimiterMatching, which accepts a string type parameterThis parameter is the file name described as above. This function repeats reading one character from the file and taking it into process until an error appears or at the end of the file. 

Submission

Submit only the following.

1. Stack.h - your completed class Stack definition 

2. Stack.cpp - if created 

3. main.cpp - your completed delimiter matching file

#include<iostream>
#include<fstream>
#include<string>
using namespace std;

template <class T>
class Stack
{
private:
	T *data;//指向栈顶指针
	int size;//栈长度
	int top;
public:
	Stack(int size=0)//构造函数
	{
		this->size=size;
		this->top=0;
		this->data=new T[size];
	}
	~Stack()//析构函数
	{
		size=0;
		top=0;
		delete []data;	
	}
	void Pop()//栈顶元素出栈
	{
		if(!isEmpty())
		{
			top--;
		}
	}
	void Push(const T & item)//入栈
	{
		if(top<size)
		{
			data[top]=item;
			top++;
		}
	}
	T & Top()//获取栈顶元素
	{
		if(!isEmpty())
		{	
			return data[top-1];
		}
	}
	bool isEmpty()
	{
		if(top==0)
		{
			cout<<"Stack is empty"<<endl;
			return true;
		}
		else return false;
	}
	int Size()
	{
		return top;
	}
};
bool delimiterMatching()
{
	cout<<"请向文件输入字符串:"<<endl;
	string s;
	cin>>s;
	ofstream out("F:\\test.txt");
	out<<s<<endl;
	out.close();//关闭文件
	ifstream get("F:\\test.txt");
	char string[50];
	get>>string;//获取文件
	int i=0;
	Stack<char> stack(100);
	while(string[i]!='\0')
	{
		//1、判断是否为左括号,是则入栈
		if(string[i]=='('||string[i]=='['||string[i]=='{')
		{
			stack.Push(string[i]);
			i++;
		}
		//2、右括号的情况
		else if(string[i]==')'||string[i]==']'||string[i]=='}')
		{
			if(stack.isEmpty())//如果栈空说明没有左匹配符
			{
				return false;
			}
			else//栈未空
			{
				char temp;
				temp=stack.Top();//获取栈顶元素,匹配的话栈顶元素出栈
				if((string[i]==')'&&temp=='(')||(string[i]=='}'&&temp=='{')||(string[i]==']'&&temp=='['))
				{
					stack.Pop();
					i++;
				}
				else return false;	
			}
		}
		//3、当读到 / 时有可能是注释符的情况
		else if(string[i]=='/')
		{
			if(string[i+1]=='/'&&stack.isEmpty())
				return true;
			else if(string[i+1]=='/'&&!stack.isEmpty())
				return false;
			else if(string[i+1]=='*')
			{
				int flag=0;
				stack.Push(string[i]);
				stack.Push(string[i+1]);
				i=i+2;        //i加2,读取/*后面的字符
				while(string[i]!='\0')
				{
					if(string[i]=='*'&&string[i+1]=='/')
					{
							stack.Pop();
							stack.Pop();
							i=i+2;
							flag=1;
							break;//中断循环
					}
					else
						i++;//未读到*,往后走
				}
				if(flag==0)
					return false;
			}
		}
		//4、读到*/说明而之前没有与之匹配的/*
		else if(string[i]=='*'&&string[i+1]=='/')
			return false;
		//5、往后遍历
		else
		{
			i++;
		}
		if(string[i]=='\0'&&stack.isEmpty())//最终栈空返回true
			return true;
		if(string[i]=='\0'&&!stack.isEmpty())//最终栈非空返回false
			return false;
	}
}

int main()
{
	Stack<int> s(10);
	/*s.Push(10);
	s.Push(20);
	cout<<s.Size()<<endl;
	s.Pop();
	cout<<s.Size()<<endl;
	cout<<s.isEmpty()<<endl;*///Test Stack

	char c='y';
	while(c=='Y'||c=='y')
	{
		if(delimiterMatching())
			cout<<"分隔符匹配成功!"<<endl;
		else
			cout<<"分隔符匹配不成功!"<<endl;
		cout<<"是否继续 (Y/N)"<<endl;
		cin>>c;
	}
	cout<<"退出程序"<<endl;
	return 0;
}


 

已知地面数据为2024年10月24日下午(具体时刻见数据文件内GPS时间),在江苏海洋大学西运动场使用SVC HR-1024i观测的无遮挡、直射遮挡标准板反射的太阳辐射辐射亮度数据(10-10 W cm-2 nm-1 sr-1),分别绘出无遮挡标准板、直射遮挡标准板反射的辐射亮度光谱曲线图。 完成过程: 3)使用上步结果,求算该时刻到达地面的太阳直射辐射被标准板反射的辐射亮度(10-10 W cm-2 nm-1 sr-1)光谱,绘出光谱曲线图。 完成过程: 4)假设标准板在所有观测波长上反射率均为1.0,求算到达地面的相对于水平面(即垂直向下)的太阳直射辐射通量密度(W m-2 nm-1)。 完成过程: 5)已知日地平均距离处太阳辐射通量密度光谱数据(fun_F0.m),数据单位是μW cm-2 nm-1,将该数据读入,将其单位转为W m-2 nm-1,绘出光谱曲线图。 完成过程: 6)使用上步结果及经验公式(日地距离校正因子可用下式计算),求算地面观测时刻,大气层外太阳辐射通量密度数据,绘出光谱曲线图。1 + 0.033 × cos( (360° × td) / 365 ) (1) 其中,td为日序数。 完成过程: 7)什么是太阳天顶角(θz),使用经验公式求算观测时刻太阳天顶角及其余弦。 cos(θz)= sin(φ) sin(δ) +cos(φ) cos(δ) cos(h)(2) 式中,φ是纬度,δ是太阳赤纬角(sin(δ)由公式3计算),h是太阳时角(由公式4计算)。 δ=-23.45摄氏度cos(360度/365(d+10)) (3) h=15(tsolar​−12) (4) 其中,d为日序数,tsolar为太阳时(由公式5近似计算)。 tsolar=tclock +4/60(λstd-λloc) (5) 其中,tclock为时区时间,λstd​是时区标准经度,λloc​是观测点经度。 完成过程: 8)使用前两步结果,求算大气层外相对于水平面(即垂直向下)的太阳辐射通量密度(W m-2 nm-1)光谱。 完成过程: 9)什么是大气质量数,使用经验公式(6)求算观测时刻大气质量数。 m = 1 / [cos(z) + 0.50572 × (96.07995° - z)^(-1.6364)](6) 完成过程: 10)使用前面步骤结果,求算大气总光学厚度光谱,绘制大气总光学厚度光谱曲线,给出550 nm处大气总光学厚度数据。 完成过程: 11)什么是瑞利散射光学厚度,已知观测时刻地面气压数据(1021.2 hPa),使用经验公式计算观测时刻大气瑞利散射光学厚度。 τ_R(λ) = (P/P₀) × (84.35/λ⁴ - 1.255/λ⁵ + 1.40/λ⁶) × 10⁻⁴(7) 其中,P为观测时刻地面气压,P0为海平面标准大气压1013.25 hPa,λ为波长(单位为微米)。 完成过程: 12)使用前面步骤结果,求算大气气溶胶光学厚度光谱,给出光谱曲线,给出550 nm AOD数据。 (可用400 ~ 680、740 ~ 753、870 ~ 883、1020 ~ 1060 nm作为大气吸收作用影响较小的波段。) 数据为gr102424_001无遮挡标准版.sig和gr102424_002遮挡标准板1.sig(位置F:\dingliang bigtest),文件结构如下: gr102424_001无遮挡标准版.sig:/** Spectra Vista SIG Data */ name= gr102424_001.sig instrument= HI: 1142031 (HR-1024i) integration= 40.0, 40.0, 10.0, 30.0, 40.0, 10.0 scan method= Time-based, Time-based scan coadds= 125, 116, 392, 166, 116, 392 scan time= 5, 5 scan settings= AI, AI external data set1= 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 external data set2= 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 external data dark= 0,0,0,0,0,0,0,0 external data mask= 0 optic= FIBER1(2), FIBER1(2) temp= 22.3, -5.0, -9.8, 22.8, -5.0, -9.8 battery= 7.51, 7.45 error= 1, 1 units= Radiance, Radiance time= 2024/10/23 9:47:11, 2024/10/23 9:48:03 longitude= 11912.7438E , 11912.7438E latitude= 3436.6458N , 3436.6453N gpstime= 061740.000 , 061832.000 comm= memory slot= 1, 3 factors= 1.200, 1.170, 1.000 [Overlap: Remove @ 990,1900, Matching Type: None] data= 336.1 19816.09 18797.35 94.86 337.7 22528.28 21318.55 94.63 339.2 24719.48 23327.15 94.37 340.7 27921.59 26424.78 94.64 342.3 30279.36 28805.83 95.13 343.8 31832.10 30374.40 95.42 345.3 31917.38 30307.43 94.96 346.8 34518.50 32761.83 94.91 … 2501.2 969.59 943.38 97.30 2503.4 1059.60 1033.11 97.50 2505.5 1079.60 1025.62 95.00 2507.6 910.80 910.80 100.00 2509.7 844.80 816.64 96.67 2511.8 772.88 744.25 96.30 2513.9 723.88 608.06 84.00 gr102424_002遮挡标准板1.sig: / Spectra Vista SIG Data ***/ name= gr102424_002.sig instrument= HI: 1142031 (HR-1024i) integration= 40.0, 40.0, 10.0, 70.0, 40.0, 10.0 scan method= Time-based, Time-based scan coadds= 125, 116, 392, 71, 116, 392 scan time= 5, 5 scan settings= AI, AI external data set1= 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 external data set2= 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 external data dark= 0,0,0,0,0,0,0,0 external data mask= 0 optic= FIBER1(2), FIBER1(2) temp= 22.3, -5.0, -9.8, 23.1, -5.0, -9.8 battery= 7.51, 7.41 error= 1, 1 units= Radiance, Radiance time= 2024/10/23 9:47:11, 2024/10/23 9:48:42 longitude= 11912.7438E , 11912.7437E latitude= 3436.6458N , 3436.6452N gpstime= 061740.000 , 061911.000 comm= memory slot= 1, 4 factors= 1.200, 1.200, 1.000 [Overlap: Remove @ 990,1900, Matching Type: None] data= 336.1 19816.09 13633.24 68.80 337.7 22528.28 15573.17 69.13 339.2 24719.48 17118.75 69.25 340.7 27921.59 19389.32 69.44 342.3 30279.36 21196.03 70.00 343.8 31832.10 22354.63 70.23 345.3 31917.38 22175.70 69.48 … 2497.0 790.97 102.06 12.90 2499.1 852.72 51.68 6.06 2501.2 969.59 52.41 5.41 2503.4 1059.60 79.47 7.50 2505.5 1079.60 134.95 12.50 2507.6 910.80 138.00 15.15 2509.7 844.80 112.64 13.33 2511.8 772.88 57.25 7.41 2513.9 723.88 202.69 28.00 (读取数据时跳过第二列) 给出每个步骤所需的Matlab函数代码以及每个函数所需的调用代码 ,fun_F0.m的结构大致如下:%https://oceancolor.gsfc.nasa.gov/docs/rsr/f0.txt %fields=wavelength,irradiance; %units=nm,uW/cm^2/nm %missing=-999; % Thuillier, G., M. Hers?, P. C. Simon, D. Labs, H. Mandel, D. Gillotay, % ! and T. Foujols, 2003, "The solar spectral irradiance from 200 % ! to 2400 nm as measured by the SOLSPEC spectrometer from the % ! ATLAS 1-2-3 and EURECA missions, Solar Physics, 214(1): 1-22 %delimiter=space; %wavelength range 200 - 2397nm at 1nm increments; function F0=fun_F0() F0=[200 0.7729 201 0.8143 202 0.8756 203 0.9263 204 1.0335 205 1.1149 206 1.1211 207 1.2612 208 1.4234 209 1.8649 210 2.6086 211 3.1691 212 3.6546 213 3.2686 214 3.9173 … 2381 6.2889 2382 6.2090 2383 6.0941 2384 5.9994 2385 5.9797 2386 5.9948 2387 6.0355 2388 6.0769 2389 6.1063 2390 6.1071 2391 6.0972 2392 6.0766 2393 6.0643 2394 6.0564 2395 6.0574 2396 6.0511 2397 6.0476]; F0(F0==-999)=NaN;
最新发布
06-30
<think>我们面对的是一个涉及多个步骤的辐射数据处理问题,需要从两个SIG文件(无遮挡和直射遮挡标准板)中读取数据,然后进行一系列计算和绘图。由于问题涉及多个步骤,我将按照步骤顺序逐一解决,并为每个步骤提供相应的Matlab函数代码和调用代码。注意:由于数据文件包含头信息和数据部分,且数据部分格式为:波长(nm)无遮挡辐射亮度遮挡辐射亮度百分比(或其他,根据文件描述,第二个文件的数据部分第三列是遮挡的辐射亮度,而第一列是波长,第二列是无遮挡的辐射亮度?但问题描述中第二个文件的数据部分跳过第二列?实际上,根据文件结构,每个文件的数据部分有四列:波长、第一组数据、第二组数据、百分比。但根据问题描述,我们需要的是辐射亮度数据。根据问题描述:无遮挡标准板文件:gr102424_001无遮挡标准版.sig遮挡标准板文件:gr102424_002遮挡标准板1.sig在无遮挡标准板文件中,数据部分有四列:波长、无遮挡辐射亮度(第一组?)、遮挡辐射亮度(第二组?)、百分比?但是文件头中units=Radiance, Radiance,说明有两组辐射亮度。然而,在问题描述中,我们只需要无遮挡标准板的辐射亮度和遮挡标准板的辐射亮度。注意,在无遮挡标准板文件中,数据部分第二列和第三列分别是什么?根据头信息,两次测量的时间不同,因此可能是两次测量的辐射亮度。但是问题要求绘出无遮挡标准板反射的辐射亮度光谱曲线图,所以应该取哪一列?仔细阅读文件头:无遮挡文件:data=后面是四列:波长,第一组辐射亮度(可能是无遮挡),第二组辐射亮度(可能是遮挡?),百分比(可能是两组之比?)。但是问题描述中要求无遮挡标准板,所以应该取第二列(第一组)?但是头信息中两次测量的时间不同,而且两次测量可能都是无遮挡?实际上,头信息中两次测量的时间分别是9:47:11和9:48:03,而且两次的经纬度几乎一样。但是,根据文件描述,无遮挡标准板文件应该只包含无遮挡的测量,而遮挡标准板文件包含遮挡的测量。再看文件头中的“memory slot”和“factors”等,可能两次测量分别是两个不同的存储位置(slot1和slot3)。但是,问题描述中并没有明确指定使用哪一次测量。然而,在遮挡标准板文件中,头信息中memoryslot=1,4,并且数据部分也有四列:波长,第一组(可能是无遮挡),第二组(遮挡),百分比。根据问题描述,无遮挡标准板文件应该使用无遮挡的测量,遮挡标准板文件应该使用遮挡的测量。但是,两个文件的数据部分都包含两组数据。我们注意到,在无遮挡标准板文件中,头信息中两次测量时间分别是061740.000和061832.000(GPS时间),而遮挡标准板文件中两次测量时间是061740.000和061911.000。其中061740.000这个时间在两个文件中都有,而且都是第一组数据。因此,我们猜测第一组数据(第一个存储位置)可能是无遮挡的参考测量(因为两个文件的第一组数据时间相同,且都是slot1),而第二组数据则是各自条件下的测量(无遮挡文件第二组是slot3,遮挡文件第二组是slot4)。但是,问题描述要求的是无遮挡标准板和遮挡标准板反射的辐射亮度。因此,在无遮挡标准板文件中,我们使用第二组数据(即无遮挡条件下的测量)?还是使用第一组数据?这里存在歧义。然而,问题描述中明确说“无遮挡标准板”和“直射遮挡标准板”,所以:无遮挡标准板文件:我们取第二组数据(因为第一组在两个文件中相同,可能是同一个参考,而第二组才是该文件对应的条件)。但是,在无遮挡标准板文件中,头信息中两次测量的单位都是Radiance,所以第二组数据(第三列)是无遮挡的辐射亮度?但是,在数据部分,第三列被描述为遮挡辐射亮度?这似乎与文件名矛盾。重新阅读数据部分示例:无遮挡文件:336.119816.0918797.3594.86 ->第三列是18797.35,第四列是94.86(百分比),这个百分比是第三列除以第二列?那么第二列可能是无遮挡,第三列是遮挡?但文件名是无遮挡标准板,为什么有遮挡的数据?同样,遮挡文件:336.119816.0913633.2468.80 ->第二列19816.09(与无遮挡文件第二列相同),第三列13633.24(遮挡),百分比68.80。因此,我们推断:两个文件的第一组数据(第二列)是相同的,都是同一个参考测量(可能是无遮挡的参考板?),而第二组数据(第三列)则是当前文件条件下的测量。所以:无遮挡标准板文件:第三列是无遮挡标准板的辐射亮度(因为百分比是94.86%,即比参考略低,可能是不同时间或不同板子的差异?)遮挡标准板文件:第三列是遮挡标准板的辐射亮度(百分比68.80%)但是,问题要求的是无遮挡标准板和直射遮挡标准板的辐射亮度,所以我们应该取两个文件的第三列作为各自条件的辐射亮度。然而,问题描述中在读取遮挡标准板文件时说“读取数据时跳过第二列”,这意味着在遮挡标准板文件中,我们只读取波长(第一列)和第三列(遮挡的辐射亮度)?那么无遮挡标准板文件呢?同样,我们只读取第一列和第三列(无遮挡的辐射亮度)?因此,我们这样读取:无遮挡文件:第一列波长,第三列辐射亮度(无遮挡)遮挡文件:第一列波长,第三列辐射亮度(遮挡)步骤3:求算该时刻到达地面的太阳直射辐射被标准板反射的辐射亮度(即直射辐射亮度)光谱。根据辐射传输原理,无遮挡标准板反射的总辐射亮度(L_total)包含直射辐射和漫射辐射,而遮挡标准板只反射漫射辐射(因为直射被遮挡)。因此,直射辐射亮度=无遮挡辐射亮度-遮挡辐射亮度。步骤4:假设标准板反射率均为1.0,求算到达地面的相对于水平面的太阳直射辐射通量密度(Wm-2 nm-1)。注意:辐射亮度(L)和辐射通量密度(E)的关系。对于朗伯体,有:E= π* L(因为朗伯体在半球空间内的积分)。所以,直射辐射通量密度 =π *直射辐射亮度(注意单位转换:原辐射亮度单位是10^{-10}W cm^{-2} nm^{-1} sr^{-1},需要转换为W m^{-2} nm^{-1},即乘以10^{-6}(因为1cm^2=10^{-4} m^2,所以乘以10^{-4},但原单位是10^{-10},所以综合为乘以10^{-6})?不对,原单位是10^{-10}W cm^{-2} nm^{-1} sr^{-1},转换为Wm^{-2}nm^{-1}sr^{-1}:1W cm^{-2} =10000 Wm^{-2},所以乘以10000,但又有10^{-10},所以乘以10^{-10}*10000=10^{-6}。然后乘以π得到通量密度(注意通量密度是单位面积上的功率,没有立体角,而朗伯体反射的辐射通量密度与辐射亮度的关系为E=πL)。因此,步骤4:E_direct =π *(L_direct) *10^{-6}[单位:W m^{-2} nm^{-1}]其中L_direct是步骤3得到的直射辐射亮度(单位:10^{-10} Wcm^{-2}nm^{-1}sr^{-1})步骤5:读取fun_F0.m中的大气层外太阳辐射通量密度(日地平均距离处),单位是μW cm^{-2} nm^{-1},转换为Wm^{-2}nm^{-1}:乘以10^{-6}(μW到W)和乘以10000(cm^{-2}到m^{-2}),即乘以10^{-2}。步骤6:使用日地距离校正因子校正得到地面观测时刻大气层外的太阳辐射通量密度。校正因子:f=1 +0.033* cos(2*pi*td /365)(注意:公式中为360度,所以用弧度制时是2*pi,但Matlab的cos函数使用弧度,所以将360°换成2*pi弧度,即:cos(2*pi *td /365))注意:公式(1)中写的是:1+0.033 ×cos((360° ×td)/365 ),所以角度是度数,因此需要转换为弧度?或者Matlab的cos函数默认是弧度,所以我们需要将度数转换为弧度?或者公式中的360°就是360度,所以直接代入度数?实际上,公式中的360°表示360度,而cos函数通常用弧度,所以应该将度数转换为弧度:cos( (360 *td /365) *pi/180),但也可以写成:cosd(360* td/365)(使用cosd函数,参数为度数)。因此,校正后的大气层外辐射通量密度:F0_obs =F0* f(注意:F0是步骤5读入并转换后的数据)步骤7:计算太阳天顶角θz及其余弦。公式(2):cos(θz) =sin(φ)sin(δ)+ cos(φ) cos(δ) cos(h)其中,φ是纬度(单位:度),δ是太阳赤纬角(单位:度),h是太阳时角(单位:度)。公式(3):δ= -23.45 *cos(360*(d+10)/365)(注意:这里360是度数,d是日序数)公式(4):h =15* (tsolar-12)(tsolar是太阳时,单位小时)公式(5):tsolar= tclock +4/60*(λstd- λloc)(tclock是时区时间,单位小时;λstd和λloc是经度,单位度,东经为正)注意:经度格式:在数据文件中,经度是11912.7438E,表示东经119度12.7438分(因为度分格式,通常度分格式是dddmm.mmmm)。所以需要转换为度:119 +12.7438/60 =119.2123967度。纬度3436.6458N表示北纬34度36.6458分,即34 +36.6458/60 =34.6107633度。时区:中国使用东八区(UTC+8),所以λstd=120°E。观测点经度λloc=119.2123967°E(从文件读取)。观测时间:根据文件头,GPS时间为061740.000(第一个测量时间),表示06时17分40.000秒,即6:17:40(UTC时间?)。注意,GPS时间与UTC时间相差整数秒(目前大约差18秒,但我们可以忽略,因为要求不高)?或者,根据问题描述,地面数据为2024年10月24日下午,但文件中的日期是2024/10/23?这可能是日期错误?需要确认。但是,我们按照文件中的日期和时间来计算。文件中的日期是2024/10/23,时间有两个:第一个测量时间(GPS时间061740.000)对应的是9:47:11(当地时间?)。这里存在时区转换问题。GPS时间是UTC时间,而文件中的time是当地时间(北京时)?因为经度是东经,所以当地时区时间(北京时间)应该是UTC+8。因此,GPS时间061740.000(UTC时间6:17:40)对应的北京时间(UTC+8)是14:17:40?但是文件中的time却是9:47:11?这似乎不一致。可能文件中的time是UTC时间?或者GPS时间就是UTC时间,而文件中的time是另一个时间?这里需要澄清。由于问题复杂,我们假设使用文件头中的gpstime(第一个测量时间)作为UTC时间,然后转换为当地时间(北京时间)?但是公式(5)中tclock是时区时间(即北京时间),所以我们需要知道观测时刻的北京时间。根据文件头中的time字段:2024/10/239:47:11(无遮挡文件的第一个测量时间)和9:48:03(第二个测量时间)。我们使用第一个测量时间(9:47:11)作为tclock(北京时间,因为江苏海洋大学在东八区)。因此,我们取tclock=9.786388889小时(9小时47分11秒=9+47/60+11/3600)。经度:λstd=120°(东八区标准经度),λloc=119.2123967°(从经度字符串转换得到)。计算太阳时:tsolar= tclock +4/60*(λstd- λloc)=9.786388889+ (4/60)*(120-119.2123967)=9.786388889+ (4/60)*0.7876033 ≈9.786388889+0.052506887≈9.838895776小时。然后计算时角:h =15*(tsolar-12)=15*(9.838895776-12)=15*(-2.161104224)=-32.41656336度。纬度φ:34.6107633度(北纬,取正值)。日序数d:2024年10月23日是一年中的第几天?注意,2024是闰年(能被4整除且不是整百年,或者能被400整除)。2024年2月有29天。计算日序数:1月31天,229天,3月31天,4月30天,5月31天,6月30天,7月31天,8月31天,9月30天,10月23日。总天数:31+29+31+30+31+30+31+31+30+23=297天(注意:1月到9月是前9个月,10月23日是第23天,所以总天数=31(1月)+29(2月)+31(3月)+30(4月)+31(5月)+30(6月)+31(7月)+31(8月)+30(9月)+23=297天)。计算太阳赤纬角δ:公式(3)δ =-23.45* cosd(360*(d+10)/365) =-23.45* cosd(360*(297+10)/365)= -23.45 *cosd(360*307/365) =-23.45* cosd(302.7945205)注意:cosd(302.7945205)=cosd(360-57.2054795)=cosd(-57.2054795)=cosd(57.2054795)≈0.5426(计算器计算),所以δ≈-23.45*0.5426≈-12.72度。然后计算cos(θz):使用公式(2),注意三角函数用度数:cos_theta_z =sind(φ)*sind(δ) +cosd(φ)*cosd(δ)*cosd(h)代入:φ=34.6107633,δ=-12.72,h=-32.41656336注意:sin和cos函数在Matlab中使用度数时用sind和cosd。步骤8:求算大气层外相对于水平面的太阳辐射通量密度(Wm^{-2}nm^{-1})光谱。大气层外太阳辐射通量密度是垂直于太阳光线方向(即法向)的,而我们要的是相对于水平面(垂直向下)的。根据几何关系,水平面上的直射辐射通量密度等于法向辐射通量密度乘以太阳天顶角的余弦(即cosθz)。所以:E_horizontal =F0_obs* cosθz(F0_obs是步骤6得到的法向大气层外辐射通量密度)步骤9:计算大气质量数m(公式6):m =1/ [cosd(z)+0.50572* (96.07995- z)^(-1.6364)](z是太阳天顶角,单位度)步骤10:求算大气总光学厚度τ_total(λ)光谱。根据比尔定律,地面水平面上的直射辐射通量密度(E_direct)与大气层外水平面上的直射辐射通量密度(E_horizontal)之间的关系为:E_direct =E_horizontal *exp(-τ_total* m)所以:τ_total(λ) =- (1/m) *ln(E_direct/E_horizontal)注意:E_direct是步骤4得到的地面水平面上的太阳直射辐射通量密度(单位Wm^{-2}nm^{-1}),E_horizontal是步骤8得到的大气层外水平面上的太阳直射辐射通量密度(单位Wm^{-2}nm^{-1}),m是步骤9得到的大气质量数(标量,与波长无关)。但是,注意:步骤4的E_direct是假设标准板反射率为1.0计算得到的,而实际上标准板反射率可能不是1.0,但我们假设为1.0,所以E_direct就是地面到达的直射辐射通量密度(在水平面上)。步骤11:计算瑞利散射光学厚度(公式7):τ_R(λ) =(P/P0) *(84.35/λ^4-1.255/λ^5 +1.40/λ^6) *10^{-4}其中,λ单位为微米(注意:我们的波长数据单位是nm,需要转换为微米:除以1000),P=1021.2hPa,P0=1013.25 hPa。步骤12:求算气溶胶光学厚度τ_a(λ):τ_total =τ_R+ τ_a +其他(如水汽吸收等),但在选定波段(400~680,740~753,870~883,1020~1060 nm)内,可以认为其他吸收作用较小,所以τ_a=τ_total- τ_R然后,我们给出550 nm处的AOD(τ_a(550))以及整个光谱的τ_a曲线。由于问题复杂,我们将每个步骤写成独立的函数,并在主程序中调用。注意:数据读取函数需要跳过文件头,读取数据部分。两个数据文件的数据部分从“data=”之后开始,每行四列(波长,第一组,第二组,百分比),我们只读取波长和第三列(第二组数据)。开始编写函数:1.读取SIG文件的函数(读取波长和辐射亮度,注意跳过第二列,即我们只取第一列和第三列):</think>以下我将逐步为您提供每个步骤的 MATLAB 实现代码和详细解释。代码严格按照问题要求编写,并包含必要的注释。 ### 数据准备函数 首先创建读取 SIG 文件和 F0 数据的函数: ```matlab function [wavelength, radiance] = read_sig_file(filename, skip_second_col) % 读取SIG文件,跳过指定列 fid = fopen(filename, 'r'); data_started = false; data_lines = {}; while ~feof(fid) line = fgetl(fid); if strcmp(strtrim(line), 'data=') data_started = true; continue; end if data_started data_lines{end+1} = line; end end fclose(fid); % 解析数据 data = zeros(length(data_lines), 4); for i = 1:length(data_lines) data(i,:) = sscanf(data_lines{i}, '%f')'; end wavelength = data(:,1); if skip_second_col radiance = data(:,3); % 跳过第二列 else radiance = data(:,2); % 不跳过第二列 end end function F0 = fun_F0() % 内置日地平均距离处太阳辐射数据 F0 = [ 200 0.7729 201 0.8143 % ... 完整数据见问题描述 2397 6.0476 ]; F0(F0(:,2) == -999, 2) = NaN; end ``` ### 主程序步骤实现 ```matlab %% 步骤1-2: 读取数据并绘制光谱曲线 [wl_unshaded, rad_unshaded] = read_sig_file('gr102424_001无遮挡标准版.sig', false); [wl_shaded, rad_shaded] = read_sig_file('gr102424_002遮挡标准板1.sig', true); % 确保波长范围一致 [common_wl, idx1, idx2] = intersect(wl_unshaded, wl_shaded); rad_unshaded = rad_unshaded(idx1); rad_shaded = rad_shaded(idx2); % 绘图 figure; plot(common_wl, rad_unshaded, 'b', 'LineWidth', 1.5); hold on; plot(common_wl, rad_shaded, 'r', 'LineWidth', 1.5); xlabel('波长 (nm)'); ylabel('辐射亮度 (10^{-10} W cm^{-2} nm^{-1} sr^{-1})'); legend('无遮挡', '直射遮挡'); title('标准板反射辐射亮度'); grid on; %% 步骤3: 计算直射辐射亮度 direct_radiance = rad_unshaded - rad_shaded; % 绘图 figure; plot(common_wl, direct_radiance, 'g', 'LineWidth', 1.5); xlabel('波长 (nm)'); ylabel('直射辐射亮度 (10^{-10} W cm^{-2} nm^{-1} sr^{-1})'); title('太阳直射辐射被标准板反射的辐射亮度'); grid on; %% 步骤4: 计算太阳直射辐射通量密度 % 单位转换: 10^{-10} W/cm&sup2; → W/m&sup2; (1 cm&sup2; = 10^{-4} m&sup2; → 乘以10^6) direct_flux_density = pi * direct_radiance * 1e-6; % 绘图 figure; plot(common_wl, direct_flux_density, 'm', 'LineWidth', 1.5); xlabel('波长 (nm)'); ylabel('辐射通量密度 (W m^{-2} nm^{-1})'); title('地面太阳直射辐射通量密度'); grid on; %% 步骤5: 读取并转换F0数据 F0_data = fun_F0(); F0_wl = F0_data(:,1); F0_irrad = F0_data(:,2) * 1e-2; % μW/cm&sup2; → W/m&sup2; (1e-6 * 10000 = 1e-2) % 绘图 figure; plot(F0_wl, F0_irrad, 'k', 'LineWidth', 1.5); xlabel('波长 (nm)'); ylabel('辐射通量密度 (W m^{-2} nm^{-1})'); title('日地平均距离处太阳辐射'); grid on; %% 步骤6: 日地距离校正 % 计算日序数 (2024年10月24日,闰年) month_days = [31, 29, 31, 30, 31, 30, 31, 31, 30]; % 1-9月 d = sum(month_days) + 24; % td = 297 (10月24日) f = 1 + 0.033 * cosd(360 * d / 365); F0_corrected = F0_irrad * f; % 绘图 figure; plot(F0_wl, F0_corrected, 'b', 'LineWidth', 1.5); xlabel('波长 (nm)'); ylabel('辐射通量密度 (W m^{-2} nm^{-1})'); title('校正后大气层外太阳辐射'); grid on; %% 步骤7: 计算太阳天顶角 % 解析经纬度 (示例: "11912.7438E" → 119 + 12.7438/60 = 119.2124°E) lon_str = '11912.7438E'; lat_str = '3436.6458N'; lon_deg = str2double(lon_str(1:3)) + str2double(lon_str(4:end-1))/60; lat_deg = str2double(lat_str(1:2)) + str2double(lat_str(3:end-1))/60; if lon_str(end) == 'W', lon_deg = -lon_deg; end if lat_str(end) == 'S', lat_deg = -lat_deg; end % 计算太阳赤纬角 (δ) delta_deg = -23.45 * cosd(360*(d + 10)/365); % 计算太阳时 (UTC+8) gpstime = '061740.000'; % HHMMSS.sss utc_hour = str2double(gpstime(1:2)) + str2double(gpstime(3:4))/60 + str2double(gpstime(5:end))/3600; t_clock = utc_hour + 8; % UTC+8 转换为本地时间 % 太阳时计算 (λ_std=120° for UTC+8) t_solar = t_clock + (4/60)*(120 - lon_deg); h_deg = 15 * (t_solar - 12); % 时角 % 计算天顶角余弦 cos_theta_z = sind(lat_deg)*sind(delta_deg) + cosd(lat_deg)*cosd(delta_deg)*cosd(h_deg); theta_z = acosd(cos_theta_z); fprintf('太阳天顶角: %.2f°\n余弦值: %.4f\n', theta_z, cos_theta_z); %% 步骤8: 大气层外水平面辐射 F0_horizontal = F0_corrected * cos_theta_z; % 绘图 figure; plot(F0_wl, F0_horizontal, 'r', 'LineWidth', 1.5); xlabel('波长 (nm)'); ylabel('辐射通量密度 (W m^{-2} nm^{-1})'); title('大气层外水平面太阳辐射'); grid on; %% 步骤9: 计算大气质量数 m = 1 / (cosd(theta_z) + 0.50572*(96.07995 - theta_z)^(-1.6364)); fprintf('大气质量数: %.4f\n', m); %% 步骤10: 计算大气总光学厚度 % 插值到相同波长 F0_horizontal_interp = interp1(F0_wl, F0_horizontal, common_wl, 'linear', 'extrap'); tau_total = -log(direct_flux_density ./ F0_horizontal_interp) / m; % 绘图 figure; plot(common_wl, tau_total, 'b', 'LineWidth', 1.5); xlabel('波长 (nm)'); ylabel('大气总光学厚度'); title('大气总光学厚度光谱'); grid on; % 550nm处结果 [~, idx550] = min(abs(common_wl - 550)); fprintf('550nm大气总光学厚度: %.4f\n', tau_total(idx550)); %% 步骤11: 瑞利散射光学厚度 P = 1021.2; % hPa P0 = 1013.25; % hPa lambda_um = common_wl / 1000; % nm → μm tau_R = (P/P0) * (84.35./lambda_um.^4 - 1.255./lambda_um.^5 + 1.40./lambda_um.^6) * 1e-4; % 绘图 figure; plot(common_wl, tau_R, 'g', 'LineWidth', 1.5); xlabel('波长 (nm)'); ylabel('瑞利光学厚度'); title('瑞利散射光学厚度'); grid on; %% 步骤12: 计算气溶胶光学厚度(AOD) tau_a = tau_total - tau_R; % 选择吸收作用较小波段 bands = [400 680; 740 753; 870 883; 1020 1060]; mask = false(size(common_wl)); for i = 1:size(bands,1) mask = mask | (common_wl >= bands(i,1) & common_wl <= bands(i,2)); end % 绘图 figure; plot(common_wl, tau_a, 'b', common_wl(mask), tau_a(mask), 'ro'); xlabel('波长 (nm)'); ylabel('气溶胶光学厚度(AOD)'); title('气溶胶光学厚度光谱'); legend('全波段', '低吸收波段', 'Location', 'best'); grid on; % 550nm处AOD fprintf('550nm气溶胶光学厚度(AOD): %.4f\n', tau_a(idx550)); ``` ### 关键步骤解释 1. **数据读取**: - `read_sig_file` 函数处理 SIG 文件格式,跳过文件头直接读取数据部分 - 单位转换包含辐射亮度(10⁻¹⁰ → W/m&sup2;)和 F0 数据(μW/cm&sup2; → W/m&sup2;) 2. **日地距离校正**: ```matlab f = 1 + 0.033 * cosd(360 * d / 365); % 校正因子 F0_corrected = F0_irrad * f; % 应用校正 ``` 3. **太阳位置计算**: - 赤纬角:`delta_deg = -23.45 * cosd(360*(d + 10)/365)` - 时角:`h_deg = 15 * (t_solar - 12)` - 天顶角:`cos_theta_z = sind(lat)*sind(delta) + cosd(lat)*cosd(delta)*cosd(h)` 4. **大气光学厚度计算**: ```matlab tau_total = -log(direct_flux_density ./ F0_horizontal) / m; tau_a = tau_total - tau_R; % 气溶胶光学厚度 ``` 5. **波段选择**: ```matlab bands = [400 680; 740 753; 870 883; 1020 1060]; % 低吸收波段 mask = (wavelength >= band_min & wavelength <= band_max); ``` ### 注意事项 1. **单位转换关键**: - 辐射亮度:10⁻¹⁰ W/cm&sup2; → 10⁻⁶ W/m&sup2; (乘以 10⁶) - F0 数据:μW/cm&sup2; = 10⁻⁶ W/cm&sup2; → W/m&sup2; 需 ×10⁻&sup2; (10⁻⁶×10⁴) 2. **坐标转换**: - 经度 "11912.7438E" → 119° + 12.7438'/60 = 119.2124° - GPS 时间需转换为本地时间(UTC+8) 3. **插值处理**: - 使用 `interp1` 将 F0 数据插值到观测波长 - 外推处理确保全波段覆盖 4. **物理模型**: - 比尔定律:`E = E0 * exp(-τm)` - 瑞利散射:τ_R ∝ λ⁻⁴ 5. **异常处理**: - 吸收波段选择避免水汽/臭氧吸收区 - NaN 数据处理(F0 中的 -999)
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值