20170803

#include <stdio.h>
#include <math.h>
#include "sofa.h"

int main ()
{
   int iy, im, id, ih, min, j;
   double sec, xp, yp, sp,dut1, ddp80, dde80, dx00, dy00, dx06, dy06,
   djmjd0, date, time, utc, dat, tai, tt, tut, ut1, rp[3][3],
   dp80, de80, dpsi, deps, epsa, rn[3][3], rnpb[3][3], ee,gst,
   rc2ti[3][3], rpom[3][3], rc2it[3][3], x, y, s, rc2i[3][3],
   era, dp00, de00, rb[3][3], rpb[3][3], v1[3], v2[3], ddp00,
   dde00,pv[2][3],pv_star[2][3],p_star[3],pvh[2][3],pvb[2][3],
   pv_earth[2][3],v_earth[3],c[3],
   elong,phi,hm,theta;
   double ra,dec, pmr,pmd,px,rv;

   iauASTROM* astrom;
   /* UTC. */
   iy = 2017;
   im = 4;
   id = 5;
   ih = 12;
   min = 0;
   sec = 0.0;
   /* Polar motion (arcsec->radians). */
   xp = 0.0349282 * DAS2R;
   yp = 0.4833163 * DAS2R;
   /* UT1-UTC (s). */
   dut1 = -0.072073685;
   /* Nutation corrections wrt IAU 1976/1980 (mas->radians). */
   ddp80 = -55.0655 * DMAS2R;
   dde80 = -6.3580 * DMAS2R;
   /* CIP offsets wrt IAU 2000A (mas->radians). */
   dx00 = 0.1725 * DMAS2R;
   dy00 = -0.2650 * DMAS2R;
   /* CIP offsets wrt IAU 2006/2000A (mas->radians). */
   dx06 = 0.1750 * DMAS2R;
   dy06 = -0.2259 * DMAS2R;
   /* TT (MJD). */
   j = iauCal2jd ( iy, im, id, &djmjd0, &date );
   if ( j < 0 ) return j;
   time = ( 60.0*(double)(60*ih + min) + sec ) / DAYSEC;
   utc = date + time;
   j= iauDat ( iy, im, id, time, &dat );
   if ( j < 0 ) return j;
   tai = utc + dat/DAYSEC;
   tt = tai + 32.184/DAYSEC;
   /* UT1. */
   tut = time + dut1/DAYSEC;
   ut1 = date + tut;
   printf("%f %fn",djmjd0,date);

   /* Observatory */
   elong=106.856666872*DD2R;
   phi=25.6529518158*DD2R;
   hm=1138.7;

   /* Star: Ra, Dec etc */
   ra=70.0*DD2R;
   dec=20.0*DD2R;
   pmr=0.0;
   pmd=0.0;
   px=1.0E-3;
   rv=0.0;
   

   theta=iauEra00(djmjd0,date);
   sp=iauSp00(djmjd0,date);
   iauPvtob ( elong, phi, hm, xp, yp, sp, theta, pv );
   iauEpv00 ( djmjd0,date,pvh,pvb);
   //iauPvstar (pvb,&ra,&dec,&pmr,&pmd,&px,&rv);
   //iauStarpv (ra,dec,pmr,pmd,px,rv,pv_star);
   //iauPv2p (pv_star,p_star);
   //iauS2xpv (1.0/DAU,DAYSEC/DAU,pv,pv_earth);
   iauS2c (ra,DPI/2.0-dec,c);
   rv=c[0]*pvb[1][0]+c[1]*pvb[1][1]+c[2]*pvb[1][2];
   rv=rv*DAU/DAYSEC/1000;
   
   
   //printf("%fn",DD2R*57.3);
   //printf("%fn",pvb[0][0]*DAU/DAYSEC/1000);// km
   //printf("%fn",pvb[0][1]*DAU/DAYSEC/1000);
   //printf("%fn",pvb[0][2]*DAU/DAYSEC/1000);
   printf("%fn",pvb[1][0]*DAU/DAYSEC/1000);//km/s
   printf("%fn",pvb[1][1]*DAU/DAYSEC/1000);
   printf("%fn",pvb[1][2]*DAU/DAYSEC/1000);
   //
   //printf("%fn",pv[0][0]/1000+pvb[0][0]*DAU/DAYSEC/1000);
   //printf("%fn",pv[0][1]/1000+pvb[0][1]*DAU/DAYSEC/1000);
   //printf("%fn",pv[0][2]/1000+pvb[0][2]*DAU/DAYSEC/1000);
   printf("%fn",pv[1][0]/1000+pvb[1][0]*DAU/DAYSEC/1000);
   printf("%fn",pv[1][1]/1000+pvb[1][1]*DAU/DAYSEC/1000);
   printf("%fn",pv[1][2]/1000+pvb[1][2]*DAU/DAYSEC/1000);

   printf("%fn",pv_star[1][0]*DAU/DAYSEC/1000);//km/s
   printf("%fn",pv_star[1][1]*DAU/DAYSEC/1000);
   printf("%fn",pv_star[1][2]*DAU/DAYSEC/1000);

   printf("%fn",p_star[0]*DAU/1000);// km
   printf("%fn",p_star[1]*DAU/1000);
   printf("%fn",p_star[2]*DAU/1000);

   printf("%fn",rv);


   return 0;
}

20150907

编译命令脚本

#！/bin/sh

gcc -c $1.c
gcc -lgsl -lgslcblas -lm -lutil $1.o
./a.out


用函数叠代定义傅里叶级数

#include <stdio.h>
#include <math.h>

double fourier(double x,int n)
{
  if(n<0)
  {
     printf("n should be equal or larger than 0");
  }
  else if(n==0)
  {
     return 1.0;
  }
  else
  {
     return fourier(x,n-1)+cos(2.0*M_PI*n*x);
  }
}


int main (void)
{
 double y;
 y=fourier(0.1,10);
 printf("% .18fn",y);
 return 0;
}