利用日期、经纬度求日出日落时间 C语言程序代码（zz）


先贴在这了，后面应该用得着
http://zhidao.baidu.com/link?url=iw-hcd_tLpRtf4r2Kh-NmDPaQ10UdlunBQUWaz14J-eNEq5fw-y83ydgOTn7KX9Z_t5rGywtlFGaygRJeEYBPa  

  1. #define PI 3.1415926 
  2.  
  3.  #include<math.h> 
  4.  
  5.  #include<iostream> 
  6.  
  7.  using namespace std; 
  8.  
  9.  
  10.  int days_of_month_1[]={31,28,31,30,31,30,31,31,30,31,30,31}; 
  11.  
  12.  int days_of_month_2[]={31,29,31,30,31,30,31,31,30,31,30,31}; 
  13.  
  14.  long double h=-0.833; 
  15.  
  16.  //定义全局变量 
  17.  
  18.  
  19. void input_date(int c[]){ 
  20.  
  21.      int i; 
  22.  
  23.      cout<<"Enter the date (form: 2009 03 10):"<<endl; 
  24.  
  25.      for(i=0;i<3;i++){ 
  26.  
  27.          cin>>c[i]; 
  28.  
  29.      } 
  30.  
  31.  } 
  32.  
  33.  //输入日期 
  34.  
  35.  
  36. void input_glat(int c[]){ 
  37.  
  38.      int i; 
  39.  
  40.      cout<<"Enter the degree of latitude(range: 0°- 60°,form: 40 40 40 (means 40°40′40″)):"<<endl; 
  41.  
  42.      for(i=0;i<3;i++){ 
  43.  
  44.          cin>>c[i]; 
  45.  
  46.      } 
  47.  
  48.  } 
  49.  
  50.  //输入纬度 
  51.  
  52.  
  53. void input_glong(int c[]){ 
  54.  
  55.      int i; 
  56.  
  57.      cout<<"Enter the degree of longitude(west is negativ,form: 40 40 40 (means 40°40′40″)):"<<endl; 
  58.  
  59.      for(i=0;i<3;i++){ 
  60.  
  61.          cin>>c[i]; 
  62.  
  63.      } 
  64.  
  65.  } 
  66.  
  67.  //输入经度 
  68.  
  69.  
  70. int leap_year(int year){ 
  71.  
  72.      if(((year%400==0) || (year%100!=0) && (year%4==0))) return 1; 
  73.  
  74.      else return 0; 
  75.  
  76.  } 
  77.  
  78.  //判断是否为闰年:若为闰年，返回1；若非闰年，返回0 
  79.  
  80.  
  81.  int days(int year, int month, int date){ 
  82.  
  83.      int i,a=0; 
  84.  
  85.      for(i=2000;i<year;i++){ 
  86.  
  87.          if(leap_year(i)) a=a+366; 
  88.  
  89.          else a=a+365; 
  90.  
  91.      } 
  92.  
  93.      if(leap_year(year)){ 
  94.  
  95.          for(i=0;i<month-1;i++){ 
  96.  
  97.              a=a+days_of_month_2[i]; 
  98.  
  99.          } 
  100.  
  101.      } 
  102.  
  103.      else { 
  104.  
  105.          for(i=0;i<month-1;i++){ 
  106.  
  107.              a=a+days_of_month_1[i]; 
  108.  
  109.          } 
  110.  
  111.      } 
  112.  
  113.      a=a+date; 
  114.  
  115.      return a; 
  116.  
  117.  } 
  118.  
  119.  //求从格林威治时间公元2000年1月1日到计算日天数days 
  120.  
  121.  
  122.  long double t_century(int days, long double UTo){ 
  123.  
  124.      return ((long double)days+UTo/360)/36525; 
  125.  
  126.  } 
  127.  
  128.  //求格林威治时间公元2000年1月1日到计算日的世纪数t 
  129.  
  130.  
  131.  long double L_sun(long double t_century){ 
  132.  
  133.      return (280.460+36000.770*t_century); 
  134.  
  135.  } 
  136.  
  137.  //求太阳的平黄径 
  138.  
  139.  
  140. long double G_sun(long double t_century){ 
  141.  
  142.      return (357.528+35999.050*t_century); 
  143.  
  144.  } 
  145.  
  146.  //求太阳的平近点角 
  147.  
  148.  
  149. long double ecliptic_longitude(long double L_sun,long double G_sun){ 
  150.  
  151.      return (L_sun+1.915*sin(G_sun*PI/180)+0.02*sin(2*G_sun*PI/180)); 
  152.  
  153.  } 
  154.  
  155.  //求黄道经度 
  156.  
  157.  
  158. long double earth_tilt(long double t_century){ 
  159.  
  160.      return (23.4393-0.0130*t_century); 
  161.  
  162.  } 
  163.  
  164.  //求地球倾角 
  165.  
  166.  
  167. long double sun_deviation(long double earth_tilt, long double ecliptic_longitude){ 
  168.  
  169.      return (180/PI*asin(sin(PI/180*earth_tilt)*sin(PI/180*ecliptic_longitude))); 
  170.  
  171.  } 
  172.  
  173.  //求太阳偏差 
  174.  
  175.  
  176. long double GHA(long double UTo, long double G_sun, long double ecliptic_longitude){ 
  177.  
  178.      return (UTo-180-1.915*sin(G_sun*PI/180)-0.02*sin(2*G_sun*PI/180)+2.466*sin(2*ecliptic_longitude*PI/180)-0.053*sin(4*ecliptic_longitude*PI/180)); 
  179.  
  180.  } 
  181.  
  182.  //求格林威治时间的太阳时间角GHA 
  183.  
  184.  
  185.  long double e(long double h, long double glat, long double sun_deviation){ 
  186.  
  187.      return 180/PI*acos((sin(h*PI/180)-sin(glat*PI/180)*sin(sun_deviation*PI/180))/(cos(glat*PI/180)*cos(sun_deviation*PI/180))); 
  188.  
  189.  } 
  190.  
  191.  //求修正值e 
  192.  
  193.  
  194.  long double UT_rise(long double UTo, long double GHA, long double glong, long double e){ 
  195.  
  196.      return (UTo-(GHA+glong+e)); 
  197.  
  198.  } 
  199.  
  200.  //求日出时间 
  201.  
  202.  
  203. long double UT_set(long double UTo, long double GHA, long double glong, long double e){ 
  204.  
  205.      return (UTo-(GHA+glong-e)); 
  206.  
  207.  } 
  208.  
  209.  //求日落时间 
  210.  
  211.  
  212. long double result_rise(long double UT, long double UTo, long double glong, long double glat, int year, int month, int date){ 
  213.  
  214.      long double d; 
  215.  
  216.      if(UT>=UTo) d=UT-UTo; 
  217.  
  218.      else d=UTo-UT; 
  219.  
  220.      if(d>=0.1) { 
  221.  
  222.          UTo=UT; 
  223.  
  224.          UT=UT_rise(UTo,GHA(UTo,G_sun(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo)))),glong,e(h,glat,sun_deviation(earth_tilt(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo)))))); 
  225.  
  226.          result_rise(UT,UTo,glong,glat,year,month,date); 
  227.  
  228.      } 
  229.  
  230.      return UT; 
  231.  
  232.  } 
  233.  
  234.  //判断并返回结果（日出） 
  235.  
  236.  
  237. long double result_set(long double UT, long double UTo, long double glong, long double glat, int year, int month, int date){ 
  238.  
  239.      long double d; 
  240.  
  241.      if(UT>=UTo) d=UT-UTo; 
  242.  
  243.      else d=UTo-UT; 
  244.  
  245.      if(d>=0.1){ 
  246.  
  247.          UTo=UT; 
  248.  
  249.          UT=UT_set(UTo,GHA(UTo,G_sun(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo)))),glong,e(h,glat,sun_deviation(earth_tilt(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo)))))); 
  250.  
  251.          result_set(UT,UTo,glong,glat,year,month,date); 
  252.  
  253.      } 
  254.  
  255.      return UT; 
  256.  
  257.  } 
  258.  
  259.  //判断并返回结果（日落） 
  260.  
  261.  
  262. int Zone(long double glong){ 
  263.  
  264.      if(glong>=0) return (int)((int)(glong/15.0)+1); 
  265.  
  266.      else return (int)((int)(glong/15.0)-1); 
  267.  
  268.  } 
  269.  
  270.  //求时区 
  271.  
  272.  
  273. void output(long double rise, long double set, long double glong){ 
  274.  
  275.      if((int)(60*(rise/15+Zone(glong)-(int)(rise/15+Zone(glong))))<10) 
  276.  
  277.          cout<<"The time at which the sun rises is "<<(int)(rise/15+Zone(glong))<<":0"<<(int)(60*(rise/15+Zone(glong)-(int)(rise/15+Zone(glong))))<<" .\n"; 
  278.  
  279.      else cout<<"The time at which the sun rises is "<<(int)(rise/15+Zone(glong))<<":"<<(int)(60*(rise/15+Zone(glong)-(int)(rise/15+Zone(glong))))<<" .\n"; 
  280.  
  281.      if((int)(60*(set/15+Zone(glong)-(int)(set/15+Zone(glong))))<10) 
  282.  
  283.          cout<<"The time at which the sun sets is "<<(int)(set/15+Zone(glong))<<": "<<(int)(60*(set/15+Zone(glong)-(int)(set/15+Zone(glong))))<<" .\n"; 
  284.  
  285.      else cout<<"The time at which the sun sets is "<<(int)(set/15+Zone(glong))<<":"<<(int)(60*(set/15+Zone(glong)-(int)(set/15+Zone(glong))))<<" .\n"; 
  286.  
  287.  } 
  288.  
  289.  //打印结果 
  290.  
  291.           
  292.  
  293. int main(){ 
  294.  
  295.      long double UTo=180.0; 
  296.  
  297.      int year,month,date; 
  298.  
  299.      long double glat,glong; 
  300.  
  301.      int c[3]; 
  302.  
  303.      input_date(c); 
  304.  
  305.      year=c[0]; 
  306.  
  307.      month=c[1]; 
  308.  
  309.      date=c[2]; 
  310.  
  311.      input_glat(c); 
  312.  
  313.      glat=c[0]+c[1]/60+c[2]/3600; 
  314.  
  315.      input_glong(c); 
  316.  
  317.      glong=c[0]+c[1]/60+c[2]/3600; 
  318.  
  319.      long double rise,set; 
  320.  
  321.      rise=result_rise(UT_rise(UTo,GHA(UTo,G_sun(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo)))),glong,e(h,glat,sun_deviation(earth_tilt(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo)))))),UTo,glong,glat,year,month,date); 
  322.  
  323.      set=result_set(UT_set(UTo,GHA(UTo,G_sun(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo)))),glong,e(h,glat,sun_deviation(earth_tilt(t_century(days(year,month,date),UTo)),ecliptic_longitude(L_sun(t_century(days(year,month,date),UTo)),G_sun(t_century(days(year,month,date),UTo)))))),UTo,glong,glat,year,month,date); 
  324.  
  325.      output(rise,set,glong); 
  326.  
  327.      system("pause"); 
  328.  
  329.      return 0; 
  330.  
  331.  }





posted @ 2016-10-09 15:30  拿走不谢  阅读(2764)  评论(0)    收藏  举报