Program sluzy do szukania minimum funkcji dopasowania w przestrzeni parametrow opisujacy orbity planet. Wykorzystuje metode gradientowa opracowana przez autora tych slow. W programie nalezy podac parametry startowe oraz nazwe pliku z danymi pomiarowymi predkosci radialnej gwiazdy, oraz plikow wyjsciowych. Mozliwe, ze zamieszcze opis samej metody, ale nie wiem dokladnie kiedy. Prgram startuje z parametrow zadanych w programie i poruszajac sie w kierunku wyznaczanym orzez tak zwany "gradient wzgledny" szuka minimum funkcji dopasowania. Minimum globalne tej funkcji oznacza najlepsze parametry orbitalne. Jednak funkcja dopasowania jest skomplikowana, wiec wynik koncowy zalezy od miejsca w przestrzeni parametrow, z ktorego wystartujemy.

---

/*GRADIENT2.c*/

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

#define DZIEN    86400     /*liczba sekund w 1 dniu*/
#define HP       0.01   /*poczatkowy 'krok gradientowy'*/
#define ALFA     1.001      /*'wspolczynnik miniecia'*/
#define MP       2         /*liczba planet*/
#define LP       11        /*liczba parametrow*/ 
#define EPSILON  0.1       /*konczaca wartosc gradientu*/
#define EPSILON2 1E-20     /* Fp-Fn */
#define N        38        /*liczba pomiarow*/
#define B        0.0001    /*dokladnosc rozwiazania r-nia Keplera*/
#define PI       3.1415926536 /*liczba pi*/

float NR(float P, float tau, float t, float e, float b, float pi);
float f(float s[], float t[], float v[], float d[], int Npom, int Lpar, float b, int Lplan, float pi);
void g(float s[],float t[],float v[],float d[], int Npom, float b, int Lplan, float grad[], float pi);

main()
{
 int i,j,numerKroku;
 float t[N],v[N],d[N],s[LP],s0[LP],gradP[LP],gradN[LP],Fmin,smin[LP],Fn,Fp,modul,modul0,h[LP];
 float V0,K[MP],P[MP],e[MP],w[MP],tau[MP];
 FILE *plik, *wynik, *gradienty;
 char *nazwa, *nazwa_w, *nazwa_g;
 
 /*wczytanie danych*/
 
 nazwa="HD114729.dat";
 plik=fopen(nazwa,"r");
 
 for (i=0; i<=N-1; ++i)
 {
  fscanf(plik,"%f%f%f",&t[i],&v[i],&d[i]);
  t[i]=t[i]*DZIEN;
 }
 
 fclose(plik);
 
 /*parametry poczatkowe*/
 
 V0=-5.5;
 K[0]=15;
 P[0]=1080*DZIEN;
 e[0]=0.19;
 w[0]=81*PI/180;
 tau[0]=539*DZIEN;
 K[1]=6;
 P[1]=10.5262*DZIEN;
 e[1]=0.016;
 w[1]=53*PI/180;
 tau[1]=101.4*DZIEN;
 
 
 /* glowna czesc programu*/
 
 s0[0]=V0;
 for (j=0; j<=MP-1; ++j)
 {
   s0[5*j+1]=K[j];
   s0[5*j+2]=P[j];
   s0[5*j+3]=e[j];
   s0[5*j+4]=w[j];
   s0[5*j+5]=tau[j];
 }
 
 for (j=0; j<=LP-1; ++j)
 {
  s[j]=s0[j];
 }

 nazwa_w="HD114729_GRADIENT_1_1.dat";
 nazwa_g="HD114729_GRADIENT_2_1.dat";
 
 wynik=fopen(nazwa_w,"w");
 gradienty=fopen(nazwa_g,"w"); 
 
 numerKroku=0;
 Fn=f(s,t,v,d,N,LP,B,MP,PI);

 V0=s[0];
   
 for (j=0; j<=MP-1; ++j)
 {
  K[j]=s[5*j+1];
  P[j]=s[5*j+2]/DZIEN;
  e[j]=s[5*j+3];
  w[j]=s[5*j+4]*180/PI;
  tau[j]=s[5*j+5]/DZIEN; 
 }

 for (j=0; j<=LP-1; ++j)
  {
   smin[j]=s[j];
  }

 Fmin=Fn;
 g(s,t,v,d,N,B,MP,gradP,PI);

 modul=0;
 for (i=0; i<=LP-1; ++i)
 {
  modul=modul+atan(s0[i]*gradP[i])*atan(s0[i]*gradP[i]);
 }
 modul=sqrt(modul);
 modul0=modul;

 for (i=0; i<=LP-1; ++i)
 {
  h[i]=HP*Fn*modul0/sqrt((s0[i]*gradP[i])*(s0[i]*gradP[i]));
 }     /* poczatkowy "krok gradientowy" */

 fprintf(wynik,"%10d%20.15f%12.5f%12.5f",numerKroku,Fn,modul,V0);
 for (j=0; j<=MP-1; ++j)
 { 
  fprintf(wynik,"%12.5f%12.5f%12.6f%12.5f%13.5f",K[j],P[j],e[j],w[j],tau[j]);
 }
 fprintf(wynik,"\n");


 fprintf(gradienty,"%10d",numerKroku);
 for (i=0; i<=LP-1; ++i)
 { 
  fprintf(gradienty,"%12.5f",atan(s0[i]*gradP[i]));
 }
 fprintf(gradienty,"\n");

 for (i=0; i<=LP-1; ++i)
 {
  s[i]=s[i]*(1-h[i]*atan(s0[i]*gradP[i])/modul);
 }
 
 Fp=Fn;
 
 Fn=f(s,t,v,d,N,LP,B,MP,PI);
 numerKroku=1;

 V0=s[0];
   
 for (j=0; j<=MP-1; ++j)
 {
  K[j]=s[5*j+1];
  P[j]=s[5*j+2]/DZIEN;
  e[j]=s[5*j+3];
  w[j]=s[5*j+4]*180/PI;
  tau[j]=s[5*j+5]/DZIEN; 
 }

 g(s,t,v,d,N,B,MP,gradN,PI);

 for (i=0; i<=LP-1; ++i)
 {
  if (gradP[i]*gradN[i]<0)
  {
   h[i]=h[i]/ALFA;
  }
 }

 for (i=0; i<=LP-1; ++i)
 {
  gradP[i]=gradN[i];
 }

 modul=0;
 for (i=0; i<=LP-1; ++i)
 {
  modul=modul+atan(s0[i]*gradP[i])*atan(s0[i]*gradP[i]);
 }
 modul=sqrt(modul);

   
 fprintf(wynik,"%10d%20.15f%12.5f%12.5f",numerKroku,Fn,modul,V0);
 for (j=0; j<=MP-1; ++j)
 { 
  fprintf(wynik,"%12.5f%12.5f%12.6f%12.5f%13.5f",K[j],P[j],e[j],w[j],tau[j]);
 }
 fprintf(wynik,"\n");


 fprintf(gradienty,"%10d",numerKroku);
 for (i=0; i<=LP-1; ++i)
 { 
  fprintf(gradienty,"%12.5f",atan(s0[i]*gradP[i]));
 }
 fprintf(gradienty,"\n");
 
 if (Fn<Fmin)
 {
  Fmin=Fn;
  for (j=0; j<=LP-1; ++j)
  {
   smin[j]=s[j];
  }
 }


 while (numerKroku<1000)
 {
  numerKroku=numerKroku+1;
  
  Fp=Fn;
 
 for (i=0; i<=LP-1; ++i)
 {
  s[i]=s[i]*(1-h[i]*atan(s0[i]*gradP[i])/modul);
 }
  Fn=f(s,t,v,d,N,LP,B,MP,PI);
  
  V0=s[0];
   
  for (j=0; j<=MP-1; ++j)
  {
   K[j]=s[5*j+1];
   P[j]=s[5*j+2]/DZIEN;
   e[j]=s[5*j+3];
   w[j]=s[5*j+4]*180/PI;
   tau[j]=s[5*j+5]/DZIEN; 
  }

  g(s,t,v,d,N,B,MP,gradN,PI);

  for (i=0; i<=LP-1; ++i)
  {
   if (gradP[i]*gradN[i]<0)
   {
    h[i]=h[i]/ALFA;
   }
  }

  for (i=0; i<=LP-1; ++i)
  {
   gradP[i]=gradN[i];
  }

  modul=0;
  for (i=0; i<=LP-1; ++i)
  {
   modul=modul+atan(s0[i]*gradP[i])*atan(s0[i]*gradP[i]);
  }
  modul=sqrt(modul);


  fprintf(wynik,"%10d%20.15f%12.5f%12.5f",numerKroku,Fn,modul,V0);
  for (j=0; j<=MP-1; ++j)
  { 
   fprintf(wynik,"%12.5f%12.5f%12.6f%12.5f%13.5f",K[j],P[j],e[j],w[j],tau[j]);
  }
  fprintf(wynik,"\n");


  fprintf(gradienty,"%10d",numerKroku);
  for (i=0; i<=LP-1; ++i)
  { 
   fprintf(gradienty,"%12.5f",atan(s0[i]*gradP[i]));
  }
  fprintf(gradienty,"\n");

  if (Fn<Fmin)
  {
   Fmin=Fn;
   for (j=0; j<=LP-1; ++j)
   {
    smin[j]=s[j];
   }
  }
 
 }

 fclose(wynik);
 fclose(gradienty);

/*----------- wyniki ------------------*/



 V0=smin[0];
   
 for (j=0; j<=MP-1; ++j)
 {
  K[j]=smin[5*j+1];
  P[j]=smin[5*j+2]/DZIEN;
  e[j]=smin[5*j+3];
  w[j]=smin[5*j+4]*180/PI;
  tau[j]=smin[5*j+5]/DZIEN; 
 }
 
 return 0;
}
/*koniec funkcji main*/

/*funkcje*/

/*funkcja obliczajaca E w czasie t metoda Newtona-Raphsona IV rzedu z dokladnoscia b*/

float NR(float P, float tau, float t, float e, float b, float pi)
{
 float MM,EE,f0,f1,f2,f3,d1,d2,d3;
 
 MM=2*pi*(t-tau)/P;
 EE=MM;
 
 f0=EE-e*sin(EE)-MM;
 f1=1-e*cos(EE);
 f2=e*sin(EE);
 f3=e*cos(EE);
 d1=-f0/f1;
 d2=-f0/(f1+d1*f2/2);
 d3=-f0/(f1+d2*f2/2+d2*d2*f3/6);
 
 while ((d3>b) || (d3<-b))
 {
  EE=EE+d3;
  f0=EE-e*sin(EE)-MM;
  f1=1-e*cos(EE);
  f2=e*sin(EE);
  f3=e*cos(EE);
  d1=-f0/f1;
  d2=-f0/(f1+d1*f2/2);
  d3=-f0/(f1+d2*f2/2+d2*d2*f3/6);
 }
 
 return EE;
}


/*funkcja obliczajaca funkcje dopasowania chi^2*/

float f(float s[], float t[], float v[], float d[], int Npom, int Lpar, float b, int Lplan, float pi)
{
 float

chi,V0,K[Lplan],P[Lplan],e[Lplan],w[Lplan],tau[Lplan],beta2,pe[Lplan],niu[Lplan],EE[Lplan];
 int i,j;
 
 chi=0;
 
 for (i=0; i<=Npom-1; ++i)
 {
  V0=s[0];
  
  for (j=0; j<=Lplan-1; ++j)
  {
   K[j]=s[5*j+1];
   P[j]=s[5*j+2];
   e[j]=s[5*j+3];
   w[j]=s[5*j+4];
   tau[j]=s[5*j+5];
  }
  
  for (j=0; j<=Lplan-1; ++j)
  {
   pe[j]=sqrt((1+e[j])/(1-e[j]));
   EE[j]=NR(P[j],tau[j],t[i],e[j],b,PI);
   niu[j]=2*atan(pe[j]*tan(EE[j]/2));
  }
  
  beta2=v[i]-V0;
  
  for (j=0; j<=Lplan-1; ++j)
  {
   beta2=beta2-K[j]*(cos(w[j]+niu[j])+e[j]*cos(w[j]));
  }
  
  beta2=(beta2/d[i])*(beta2/d[i]);
  chi=chi+beta2;
 }
 
 chi=chi/(Npom-Lpar-1);
 
 return chi;
} 
   
 
/*funkcja liczaca gradient funkcji chi^2 w przestrzeni parametrow s*/

void g(float s[],float t[],float v[],float d[], int Npom, float b, int Lplan, float grad[], float pi)
{
 float V0,K[Lplan],P[Lplan],e[Lplan],w[Lplan],tau[Lplan],beta2,beta;
 float pe[Lplan],niu[Lplan],EE[Lplan],dniudP[Lplan],dniudE[Lplan],dniudtau[Lplan],lp;
 int i,j;

 lp=5*Lplan+1;

 for (j=0; j<=lp-1; ++j)
 {
  grad[j]=0;
 }
 
 for (i=0; i<=Npom-1; ++i)
 {
  V0=s[0];
  
  for (j=0; j<=Lplan-1; ++j)
  {
   K[j]=s[5*j+1];
   P[j]=s[5*j+2];
   e[j]=s[5*j+3];
   w[j]=s[5*j+4];
   tau[j]=s[5*j+5];
  }
  
  for (j=0; j<=Lplan-1; ++j)
  {
   pe[j]=sqrt((1+e[j])/(1-e[j]));
   EE[j]=NR(P[j],tau[j],t[i],e[j],b,PI);
   niu[j]=2*atan(pe[j]*tan(EE[j]/2));
   dniudE[j]=pe[j]/((cos(EE[j]/2))*(cos(EE[j]/2))+pe[j]*pe[j]*(sin(EE[j]/2))*(sin(EE[j]/2)));
   dniudP[j]=-dniudE[j]*2*pi*(t[i]-tau[j])/(P[j]*P[j]*(1-e[j]*cos(EE[j])));
   dniudtau[j]=-dniudE[j]*2*pi/(P[j]*(1-e[j]*cos(EE[j])));
  }
  
  beta2=v[i]-V0;

  for (j=0; j<=Lplan-1; ++j)
  {
   beta2=beta2-K[j]*(cos(w[j]+niu[j])+e[j]*cos(w[j]));
  }
  beta=2*beta2/(d[i]*d[i]);

  grad[0]=grad[0]-beta;                                                           /* dchi/dV0 */
 
  for (j=0; j<=Lplan-1; ++j)
  {
   grad[5*j+1]=grad[5*j+1]-beta*(cos(w[j]+niu[j])+e[j]*cos(w[j]));                /* dchi/dK(j) */
   grad[5*j+2]=grad[5*j+2]+beta*K[j]*sin(w[j]+niu[j])*dniudP[j];                  /* dchi/dP(j) */
   grad[5*j+3]=grad[5*j+3]-beta*K[j]*cos(w[j]);                                   /* dchi/de(j) */
   grad[5*j+4]=grad[5*j+4]+beta*K[j]*(sin(w[j]+niu[j])+e[j]*sin(w[j]));           /* dchi/dw(j) */
   grad[5*j+5]=grad[5*j+5]+beta*K[j]*sin(w[j]+niu[j])*dniudtau[j];                /* dchi/dtau(j) */
  }
 }  

}

-- CezaryMigaszewski - 26 Feb 2005

---

Daje mnie Naruszenie ochrony pamiêci :(.

Mozna dac inny warynek konczacy petle, np. okreslona liczba krokow (juz zrobione). Kod zostanie jeszcze poprawiony w przyszlosci (niedalekiej)

Brakuje plik wejść HD114729.dat - Też mogłybyś zrobić test żeby uważać użytkownik jeśli plik nie istnieje - jest to więcej ,,user-friendly" niż błęd Naruszenie ochrony pamięci np:

 nazwa="HD114729.dat";
 plik=fopen(nazwa,"r");

  if(dane==NULL)
    {
      printf("Niestety, nie widze plik  HD114729.dat . Prosze go przygotowac.\n");
      exit (0);
    };

Zrobiłem soft link ln -s HD114729_modelowa.dat HD114729.dat z plik od SymulowaneDane i teraz robi coś...