sun.c

每天太阳天顶角度和太阳方位角度的计算

#include <stdio.h>
#include <math.h>
#define PI 3.1415926

void Usage();
int DayNumber(int year, int month, int day);
float SunLatitude(int nDayNumber);
float TimeDifference( int nDayNumber);
float TimeNumber(int hour, int minute, int second);
float SunHeightAngle(float fLatitude, float fSunLatitude, float fTimeNumber);
float SunPosition(float fLatitude, float fSunHeightAngle, float fSunLatitude);
int FileDeal( char *cFileInput, float fLongtitude, float fLatitude, FILE *fpFileOutput);
void Usage()
{
   printf("\nUsage: sun <FileName> <OutputFileName>");
   printf("\nThe Content of <FileName>:\n<Site longitude>\n<Site latitude>\n");
   printf("year month date hour minute second\n");
   printf("year month date hour minute second\n");
   printf("\n                 ......          \n");

}


/* Calculate the daynumber in a year, for example, 04/17/97
   Day Number is 107    					*/

int DayNumber(int year, int month, int day)
{
   int nMonthDay[12] = {31,28,31,30,31,30,31,31,30,31,30,31};
   int nDayNumber = 0,count;

   if (year>60) year= 1900 + year;
   else  year = 2000 + year;

   if(year%4==0) nMonthDay[1] = 29;


   for (count = 0; count < month-1; count ++)
   {
      nDayNumber =nDayNumber + nMonthDay[count];
   }
   return nDayNumber+day;
}


/* Calculate the Sun Latitude according the Day Number calculated
   by function DayNumber                                          */

float SunLatitude(int nDayNumber)
{
float fSunAngle;
float fSunLatitude;

fSunAngle=2.*PI*(nDayNumber-1.)/365.0;

fSunLatitude=0.006894 - 0.399512*cos(fSunAngle) + 0.072075*sin(fSunAngle)
	- 0.006799*cos(2.*fSunAngle) + 0.000896*sin(2.*fSunAngle)
	- 0.002689*cos(3.*fSunAngle) + 0.001516*sin(3.*fSunAngle);

fSunLatitude = fSunLatitude*180./PI;
return fSunLatitude;
}

/* Caculate the TimeDifference from Day Number, the unit is rad.,
   fTimeDifference*229.183 is to convert rad. to minute           */

float TimeDifference( int nDayNumber)
{
float fSunAngle;
float fTimeDifference;

fSunAngle=2.*PI*(nDayNumber-1.)/365.0;

fTimeDifference = 0.000043 + 0.002061*cos(fSunAngle) - 0.032040*sin(fSunAngle)
		  - 0.014974*cos(2.*fSunAngle) - 0.040685*sin(2.*fSunAngle);
fTimeDifference  = fTimeDifference * 229.183;
return fTimeDifference;
}



float TimeNumber(int hour, int minute, int second)
{
   float fTimeNumber;
   fTimeNumber = hour + minute/60.0 + second/3600.0;
   return fTimeNumber;
}

float SunHeightAngle(float fLatitude, float fSunLatitude, float fTimeNumber)
{
   float fTimeAngle;
   float fSunHeightAngle;
   float fSinSunHeightAngle;

   if (fTimeNumber <= 12)
      fTimeAngle = (12.0 - fTimeNumber) * 15.0;
   else
      fTimeAngle = (fTimeNumber - 12.0) * 15.0;

   fTimeAngle = fTimeAngle * PI/180;

   fSinSunHeightAngle = sin(PI*fLatitude/180)*sin(PI*fSunLatitude/180) +
       cos(PI*fLatitude/180)*cos(PI*fSunLatitude/180)*cos(fTimeAngle);
   fSunHeightAngle = asin(fSinSunHeightAngle)*180.0/PI;

   return fSunHeightAngle;
}

float SunPosition(float fLatitude, float fSunHeightAngle, float fSunLatitude)
{
   float fSunPosition;
   float fCosSunPosition;

   fCosSunPosition = tan(PI*fSunHeightAngle/180.0)*tan(PI*fLatitude/180.0) -
       sin(PI*fSunLatitude/180.0)/((cos(PI*fSunHeightAngle/180))*(cos(PI*fLatitude/180.0)));
   fSunPosition = acos(fCosSunPosition);
   fSunPosition = fSunPosition * 180/PI;
   return fSunPosition;
}

int FileDeal( FILE *fpFileInput, float fLongtitude,float fLatitude, FILE *fpFileOutput)
{

float fSunLatitude1,fSunLatitude2,fSunLatitude;
float fTimeNumber, fTimeNumber1,fSunHeightAngle, fSunPosition, fTimeDifference;
int nDayNumber;
int year,month,day;
int hour,minute,second;
int i,j;


while (!feof (fpFileInput))
{
	fscanf(fpFileInput,"%d",&year);
	fscanf(fpFileInput,"%d",&month);
	fscanf(fpFileInput,"%d",&day);
	fscanf(fpFileInput,"%d",&hour);
	fscanf(fpFileInput,"%d",&minute);
	fscanf(fpFileInput,"%d",&second);
	if(feof(fpFileInput)) break;

	nDayNumber = DayNumber(year,month,day);
	fTimeNumber = TimeNumber(hour,minute,second);

	fSunLatitude1 = SunLatitude(nDayNumber);
	fSunLatitude2 = SunLatitude(nDayNumber + 1);
	fSunLatitude = (fSunLatitude2 + fSunLatitude1)/2.0;

	fTimeDifference = TimeDifference(nDayNumber);
	fTimeNumber1 = fTimeNumber + fTimeDifference/60 + (fLongtitude - 120)*4/60;
	fSunHeightAngle = SunHeightAngle(fLatitude, fSunLatitude, fTimeNumber1);
	fSunPosition = SunPosition( fLatitude, fSunHeightAngle, fSunLatitude);

	if(fTimeNumber1>12)
		fSunPosition = 180+fSunPosition;
	else
		fSunPosition = 180-fSunPosition;

	fprintf(fpFileOutput, "%10.2f %02d/%02d/%02d %02d:%02d:%02d %6.2f%6.2f %9.4f %9.4f \n",fTimeNumber1,
		year,month,day,hour,minute,second,fLongtitude,fLatitude,90-fSunHeightAngle,
		fSunPosition);
}
return 1;
}

int main( int argc, char *argv[])
{
FILE *fpFileName,*fpFileOutput;
char cInputFile[60];
float fLatitude,fLongtitude;
if (argc!=3)
   {
      Usage();
      exit(0);
   }

if ((fpFileName=fopen(argv[1],"rt"))==NULL)
{
   printf("\n File %s Cannot Be Found",argv[1]);
   exit(0);
}

if((fpFileOutput=fopen(argv[2],"wt"))==NULL)
{
   printf("\n File %s Cannot Be Created",argv[2]);
   exit(0);
}

   fscanf(fpFileName,"%f",&fLongtitude);
   fscanf(fpFileName, "%f",&fLatitude);

   FileDeal(fpFileName,fLongtitude,fLatitude,fpFileOutput);

fclose(fpFileName);
fclose(fpFileOutput);
}