QQ登录

只需一步,快速开始

太阳位置计算c程序

[复制链接]
hubble 发表于 2004-10-8 21:24 | 显示全部楼层 |阅读模式 来自: 中国–吉林–四平 联通

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?加入牧夫(请注明天文爱好者,否则无法通过审核,请勿使用gmail/outlook/aol/icloud邮箱注册)

×
根据Paul Schlyter, Stockholm, Sweden中等精度的行星位置计算方法,实现了其中太阳位置的计算方法。

/////////////////////////////////////////////////////////////////////////////////////////////////////
//                           名称:太阳位置计算程序
//                           作者:胡铂(http://hubble.lamost.org)
//                           日期  2004-09-30
//                           说明:根据Paul Schlyter, Stockholm, Sweden的中等精度计算方法实现
////////////////////////////////////////////////////////////////////////////////////////////////////
#include "stdio.h"
#include "math.h"
#define Pi  3.14159265
#define DE Pi/180
/////////////////////////////////////////////////////////////////
float jde(int Y,int M,int D,int hour,int min,int sec)
{
int f,g;
double mid1,mid2,J,JDE,A;
 if(M>=3)
 {
  f=Y;
  g=M;
 }
 if(M==1||M==2)
 {
  f=Y-1;
  g=M+12;
 };
 mid1=floor(365.25*f);
 mid2=floor(30.6001*(g+1));
 A=2-floor(f/100)+floor(f/400);
 J=mid1+mid2+D+A+1720994.5;
 JDE=J+hour/24+min/1440+sec/86400;
 return JDE;
}

/////////////////////////////////////变量声明///////////////////////////////
void main()
{ int i,year,month,day,hour,min,sec;
 double d,w,a,e,M,oblecl,L,E,xe,ye,r,v,lon,x,y,z,xequt,yequt,zequt,dist,RA,Decl;
// scanf("%d,%d,%d,%d,%d,%d",&year,&month,&day,&hour,&min,&sec);
// printf("\n");
///////////////////////////////测试数据////////////////////////////////////////////////////
year=1990;
month=4;
day=19;
hour=0;
min=0;
sec=0;
////////////////////////////////////////////轨道根数//////////////////////////////////////
 d=jde(year,month,day,hour,min,sec)-2451543.5;//相对儒略日;
  w=282.9404+4.70935*0.00001*d;//升交点经度
   a=1;
    e=0.016709-1.151*0.000000001*d;//偏心率
    M = 356.0470 + 0.9856002585 * d;//平近点角
     oblecl = 23.4393-3.563*0.0000001 * d;//黄赤交角
////////////////////////////////////////////////////////////////////////////////////////////      
 L=w+M;//太阳的平均精度
 L=fmod(L,360);
 E=M + (180/Pi) * e * sin(M*DE) * (1 + e * cos(M*DE));//(开普勒方程的近似解)
 E=fmod(E,360);
 xe=cos(E*DE) - e;//椭圆轨道上的直角坐标x
 ye=sin(E*DE) * sqrt(1 - e*e);//椭圆轨道上的直角坐标y
 r=sqrt(xe*xe+ye*ye);//距离
 v=atan2(ye,xe)*180/Pi;//真近点角
 //lon=fmod((v+w),360);//太阳的精度
 //太阳的黄道直角坐标
 lon=v+w;
 x=r*cos(lon*Pi/180);
 y=r*sin(lon*Pi/180);
 z=0;
 //太阳的赤道直角坐标
 xequt=x;
 yequt=y*cos(oblecl*Pi/180);
 zequt = y* sin(oblecl*Pi/180);
 //日地距离、赤经赤纬
 dist=sqrt(xequt*xequt+yequt*yequt);
 RA=atan2(yequt,xequt)*180/Pi;
 RA=fmod(RA,360);
 Decl=asin(zequt/r)*180/Pi;
 Decl=fmod(Decl,360);
 ////////////////////////////////////////////////////////////////////
 printf("d=%f\n",d);  
  printf("w=%f\n",w);
   printf("e=%f\n",e);
    printf("M=%f\n",M);
     printf("L=%f\n",L);
      printf("E=%f\n",E);
      printf("oblecl=%f",oblecl);
       printf("xe=%f,ye=%f\n",xe,ye);
        printf("r=%f,v=%f\n",r,v);
  printf("lon=%f\n",lon);
   printf("x=%f,y=%f\n",x,y);
    printf("xequt=%f,yequt=%f,zequt=%f\n",xequt,yequt,zequt);
     printf("dist=%f\n",dist);
      printf("RA=%f,Decl=%f\n",RA,Decl);



 


   
 
 

}
nngs 发表于 2004-10-9 01:45 | 显示全部楼层 来自: 美国 乔治亚州立大学
请问hubble, 你是否打算发一个系列的计算程序?
如果那样,建议你把这些计算程序都跟贴在一个主贴里。方便大家检索查询
回复 顶~ 砸~

使用道具 举报

 楼主| hubble 发表于 2004-10-9 09:37 | 显示全部楼层 来自: 中国–吉林–四平 联通

好主意

我确实打算写一个系列,但是不一定什么时候完成,下次发到一个帖子里面。另外已经发的是否可以挪到一个帖子里面去?
回复 顶~ 砸~

使用道具 举报

nngs 发表于 2004-10-9 10:15 | 显示全部楼层 来自: 美国–新泽西州–伯灵顿 Comcast有线通信控股股份有限公司
嗯,具体操作是要想想.

我建议你可以以这个贴为主贴,请大家不要跟贴.记下这主贴连接的地址.以后你再有新程序时,你可以新发一贴,在贴中注明这个主贴的连接.然后再在主贴中加入新贴的连接.

这样,网友不论读到这个系列的哪个贴子,都可以很容易找到主贴,进而找到其他分贴.

抱歉,这个方法麻烦了些,但即保证了论坛的特点(网友可以自由在这个系列的分贴里探讨),又容易检索.

主贴我给你加精华
回复 顶~ 砸~

使用道具 举报

 楼主| hubble 发表于 2004-10-9 16:48 | 显示全部楼层 来自: 中国–吉林–四平 联通
我已经把发的帖子都贴到了那个月球计算的帖子后面了。
回复 顶~ 砸~

使用道具 举报

福建ltycho 发表于 2004-11-1 19:48 | 显示全部楼层 来自: 中国–北京–北京 中国传媒大学
这个程序其实不难,有计算的公式就差不多了。
不过开发成图形界面的回更好:)
回复 顶~ 砸~

使用道具 举报

macross11 发表于 2004-11-1 23:30 | 显示全部楼层 来自: 中国–四川–成都 电信
的确,强烈要求图形界面。
回复 顶~ 砸~

使用道具 举报

heihoo 发表于 2004-11-4 16:18 | 显示全部楼层 来自: 中国–江苏–南京–玄武区 电信
能不能,贴个具体计算,月亮的合朔,以及24节气的计算机的例子,也好叫我们这些门外汉,轻松入门。
回复 顶~ 砸~

使用道具 举报

霍曼 发表于 2010-7-31 01:14 | 显示全部楼层 来自: 中国–北京–北京 歌华有线
太阳位置计算程序(2006-12-24 15:07:26)
转载   分类:格物偶得
/////////////////////////////////////////////////////////////////////////////////////////////////////
// 名称:太阳位置计算程序
// 作者:胡铂(http://leohubble.go.3322.org
// 日期 2004-09-30
// 说明:根据Paul Schlyter, Stockholm, Sweden的中等精度计算方法实现
////////////////////////////////////////////////////////////////////////////////////////////////////
#include "stdio.h"
#include "math.h"
#define Pi 3.14159265
#define DE Pi/180
/////////////////////////////////////////////////////////////////
float jde(int Y,int M,int D,int hour,int min,int sec)
{
int f,g;
double mid1,mid2,J,JDE,A;
if(M>=3)
{
f=Y;
g=M;
}
if(M==1||M==2)
{
f=Y-1;
g=M+12;
};
mid1=floor(365.25*f);
mid2=floor(30.6001*(g+1));
A=2-floor(f/100)+floor(f/400);
J=mid1+mid2+D+A+1720994.5;
JDE=J+hour/24+min/1440+sec/86400;
return JDE;
}
/////////////////////////////////////变量声明///////////////////////////////
void main()
{ int i,year,month,day,hour,min,sec;
double d,w,a,e,M,oblecl,L,E,xe,ye,r,v,lon,x,y,z,xequt,yequt,zequt,dist,RA,Decl;
// scanf("%d,%d,%d,%d,%d,%d",&year,&month,&day,&hour,&min,&sec);
// printf("\n");
///////////////////////////////测试数据////////////////////////////////////////////////////
year=1990;
month=4;
day=19;
hour=0;
min=0;
sec=0;
////////////////////////////////////////////轨道根数//////////////////////////////////////
d=jde(year,month,day,hour,min,sec)-2451543.5;//相对儒略日;
w=282.9404+4.70935*0.00001*d;//升交点经度
a=1;
e=0.016709-1.151*0.000000001*d;//偏心率
M = 356.0470 + 0.9856002585 * d;//平近点角
oblecl = 23.4393-3.563*0.0000001 * d;//黄赤交角
////////////////////////////////////////////////////////////////////////////////////////////
L=w+M;//太阳的平均精度
L=fmod(L,360);
E=M + (180/Pi) * e * sin(M*DE) * (1 + e * cos(M*DE));//(开普勒方程的近似解)
E=fmod(E,360);
xe=cos(E*DE) - e;//椭圆轨道上的直角坐标x
ye=sin(E*DE) * sqrt(1 - e*e);//椭圆轨道上的直角坐标y
r=sqrt(xe*xe+ye*ye);//距离
v=atan2(ye,xe)*180/Pi;//真近点角
//lon=fmod((v+w),360);//太阳的精度
//太阳的黄道直角坐标
lon=v+w;
x=r*cos(lon*Pi/180);
y=r*sin(lon*Pi/180);
z=0;
//太阳的赤道直角坐标
xequt=x;
yequt=y*cos(oblecl*Pi/180);
zequt = y* sin(oblecl*Pi/180);
//日地距离、赤经赤纬
dist=sqrt(xequt*xequt+yequt*yequt);
RA=atan2(yequt,xequt)*180/Pi;
RA=fmod(RA,360);
Decl=asin(zequt/r)*180/Pi;
Decl=fmod(Decl,360);
////////////////////////////////////////////////////////////////////
printf("d=%f\n",d);
printf("w=%f\n",w);
printf("e=%f\n",e);
printf("M=%f\n",M);
printf("L=%f\n",L);
printf("E=%f\n",E);
printf("oblecl=%f",oblecl);
printf("xe=%f,ye=%f\n",xe,ye);
printf("r=%f,v=%f\n",r,v);
printf("lon=%f\n",lon);
printf("x=%f,y=%f\n",x,y);
printf("xequt=%f,yequt=%f,zequt=%f\n",xequt,yequt,zequt);
printf("dist=%f\n",dist);
printf("RA=%f,Decl=%f\n",RA,Decl);
}
请教一下 我想得到wgs-84的坐标 怎么转换》
回复 顶~ 砸~

使用道具 举报

本版积分规则

APP下載|手机版|爱牧夫天文淘宝店|牧夫天文网 ( 公安备案号21021102000967 )|网站地图|辽ICP备19018387号

GMT+8, 2025-4-8 00:43 , Processed in 0.098334 second(s), 18 queries , Gzip On, Redis On.

Powered by Discuz! X3.5 Licensed

Copyright © 2001-2020, Tencent Cloud.

快速回复 返回顶部 返回列表