sunzenithpre.pro

根据输入的经度、纬度以及年、月、日、时间

;是否是润年
pro IsAdditionalYear, inYear, otAddYear
	otAddYear = 0	;1润年，0平年
	IF (inYear Mod 400 EQ 0) || ((inYear Mod 4 EQ 0) && (inYear Mod 100 NE 0)) THEN otAddYear = 1
end

;计算任意一天的积日
pro AccumulateDayNumber, inYear, inMonth, inDay, otAccuNum
	otAccuNum = 0
	;计算前所有月的总天数
	FOR  ii = 1, inMonth - 1 DO BEGIN
		IF ii NE 2 THEN BEGIN			;非2月要么是30天要么是31天
			otAccuNum = otAccuNum + 30 + (inMonth Mod 2)
		ENDIF ELSE BEGIN				;是2月平年28天润年29天
			IsAdditionalYear, inYear, otAddYear
			IF otAddYear EQ 1 THEN otAccuNum = otAccuNum + 29
			IF otAddYear EQ 0 THEN otAccuNum = otAccuNum + 28
		ENDELSE
	ENDFOR

	;加上这个月的今天
	otAccuNum = otAccuNum + inDay
end

pro sunposition, year,month,day,hour,minute,second,latitude,longitude, otSun

	;%daynum=daynumber(year,month,day);%计算天数
	;%sunlat=sunlatitude(year,month,day); %计算太阳纬度
	;%timediff=timediffererce(year,month,day);%计算时差
	;%timenum=timenumber(hour,minute,second);%计算时间
	;%sunheight=sunheightangle(latitude,hour,minute,second);%计算太阳高度角/方位角
	;%sunposi=sunpositon(latitude,sunheight,hour,minute,second,year,month,day);%计算太阳方位角
	;%result=filedeal(year,month,day,hour,minute,second,latitude,longitude); %计算太阳角

	AccumulateDayNumber, year, month, day, daynum

	sunangle1=2*!pi*(daynum-1)/365.;
	sunangle2=2*!pi*(daynum)/365.;
	sunlat1=0.006894 - 0.399512*cos(sunangle1) + 0.072075*sin(sunangle1)- 0.006799*cos(2.*sunangle1) + 0.000896*sin(2.*sunangle1)- 0.002689*cos(3.*sunangle1) + 0.001516*sin(3.*sunangle1);
	sunlat2=0.006894 - 0.399512*cos(sunangle2) + 0.072075*sin(sunangle2)- 0.006799*cos(2.*sunangle2) + 0.000896*sin(2.*sunangle2)- 0.002689*cos(3.*sunangle2) + 0.001516*sin(3.*sunangle2);
	sunlat=(sunlat1+sunlat2)/2.;
	sunlat=sunlat*180./!pi;

	sunangle=2*!pi*(daynum-1)/365;
	timediff1= 0.000043 + 0.002061*cos(sunangle) - 0.032040*sin(sunangle)- 0.014974*cos(2.*sunangle) - 0.040685*sin(2.*sunangle);
	timediff = timediff1 * 229.183;


	timenum1 = hour + minute/60. + second/3600.;
	timenum=timenum1+timediff/60.+(longitude-120.)*4./60.;

	 if timenum le 12. then begin
	      timeang1 = (12. - timenum) * 15;
	 endif else begin
	      timeang1 = (timenum - 12) * 15;
	 endelse
	  timeang = timeang1*!pi/180.;
	  sunheightangle1 = asin(sin(!pi*latitude/180.)*sin(!pi*sunlat/180.)+cos(!pi*latitude/180.)*cos(!pi*sunlat/180.)*cos(timeang));
	  sunheightangle= sunheightangle1*180/!pi;
	  sunzenith=90-sunheightangle;

	sunposi1=acos(tan(sunheightangle1)*tan(!pi*latitude/180.)-sin(!pi*sunlat/180.)/((cos(sunheightangle1))*(cos(!pi*latitude/180.))));
	sunposi = sunposi1*180./!pi;
	if(timenum gt 12) then begin
	     sunposi = 180.+sunposi;
	endif else begin
	     sunposi = 180.-sunposi;
	endelse

	otSun=[0.,0.];
	otSun[0]=sunzenith;
	otSun[1]=sunposi;
end

pro test
	year = 2008
	month = 6
	day = 23
	;hour = 1
	minute = 0
	second = 0
	latitude = 65
	longitude = 125

	FOR ii = 1, 24 DO BEGIN
		sunposition, year,month,day,ii,minute,second,latitude,longitude, otSun
		Print, '**', ii,': ', otSun
	ENDFOR
end

pro CalculateSun, inTLB
	strZenith = '-' & strAzimuth = '-'
	GetValueFromWidget, inTLB, otLong, otLati, otYear, otMonth, otDay, otHour, otMinute, otSecond, otValid
	IF otValid EQ 1 THEN BEGIN
		sunposition, otYear, otMonth, otDay, otHour, otMinute, otSecond, otLati, otLong, otSun
		IF N_ELEMENTS(otSun) EQ 2 THEN BEGIN
			strZenith = String(StrTrim(otSun[0], 2)) & strAzimuth = String(StrTrim(otSun[1], 2))
		ENDIF
	ENDIF

	id = WIDGET_INFO(inTLB, FIND_BY_UNAME = 'WID_TEXT_ZENITH') & IF id EQ 0 THEN RETURN
	WIDGET_CONTROL, id, SET_VALUE = strZenith

	id = WIDGET_INFO(inTLB, FIND_BY_UNAME = 'WID_TEXT_AZIMUTH') & IF id EQ 0 THEN RETURN
	WIDGET_CONTROL, id, SET_VALUE = strAzimuth
end

pro GetValueFromWidget, inTLB, otLong, otLati, otYear, otMonth, otDay, otHour, otMinute, otSecond, otValid
	otValid = 0

	id = WIDGET_INFO(inTLB, FIND_BY_UNAME = 'WID_TEXT_LONG')
	WIDGET_CONTROL, id, GET_VALUE = otLong & IF StrLen(otLong) EQ 0 THEN RETURN & IF StrTrim(otLong, 2) EQ '-' THEN RETURN
	otLong = Float(otLong)
	IF (otLong GT 180) || (otLong LT -180) THEN RETURN

	id = WIDGET_INFO(inTLB, FIND_BY_UNAME = 'WID_TEXT_LATI')
	WIDGET_CONTROL, id, GET_VALUE = otLati & IF StrLen(otLati) EQ 0 THEN RETURN & IF StrTrim(otLati, 2) EQ '-' THEN RETURN
	otLati = Float(otLati)
	IF (otLati GT 90) || (otLati LT -90) THEN RETURN

	id = WIDGET_INFO(inTLB, FIND_BY_UNAME = 'WID_TEXT_YEAR')
	WIDGET_CONTROL, id, GET_VALUE = otYear & IF StrLen(StrTrim(otYear, 2)) NE 4 THEN RETURN & IF StrTrim(otYear, 2) EQ '-' THEN RETURN
	otYear = Fix(otYear)
	IF otYear LE 0 THEN RETURN

	id = WIDGET_INFO(inTLB, FIND_BY_UNAME = 'WID_DROPLIST_MONTH')
	otMonth = WIDGET_INFO(id, /DROPLIST_SELECT) & otMonth = otMonth + 1

	id = WIDGET_INFO(inTLB, FIND_BY_UNAME = 'WID_DROPLIST_DAY')
	otDay = WIDGET_INFO(id, /DROPLIST_SELECT) & otDay = otDay + 1

	id = WIDGET_INFO(inTLB, FIND_BY_UNAME = 'WID_DROPLIST_HOUR')
	otHour = WIDGET_INFO(id, /DROPLIST_SELECT) & otHour = otHour

	id = WIDGET_INFO(inTLB, FIND_BY_UNAME = 'WID_DROPLIST_MINUTE')
	otMinute = WIDGET_INFO(id, /DROPLIST_SELECT) & otMinute = otMinute

	id = WIDGET_INFO(inTLB, FIND_BY_UNAME = 'WID_DROPLIST_SECOND')
	otSecond = WIDGET_INFO(id, /DROPLIST_SELECT) & otSecond = otSecond


	;id = WIDGET_INFO(inTLB, FIND_BY_UNAME = 'WID_TEXT_ZENITH')
	;id = WIDGET_INFO(inTLB, FIND_BY_UNAME = 'WID_TEXT_AZIMUTH')


	otValid = 1
end