Program sluzy do rysowania rzutu funkcji dopasowania na poszczegolne parametry. Funkcja dopasowania zalezy od pewnej liczby parametrow, amplitudy predkosci radialnych bedace wynikiem poszczegolnych planet, ich okresow obiegu, ekscentrycznosci, dlugosci pericentrum, czasu przejscia przez pericentrum oraz predkosci srodka masy ukladu. Funkcja dopasowania mowi nam jak dobrze dopasowana jest krzywa modelowa predkosci radialnej (ktorej parametrami sa wymienione powyzej) do danych pomiarowy. Im wartosc funkcji jest mniejsza tym lepsze mamy dopasowanie. W odpowiednich miejscach kodu nalezy podac parametry oraz nazwe pliku z danymi oraz plikow wyjsciowych.



/*chi2.c*/

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

#define DZIEN    86400     /*liczba sekund w 1 dniu*/
#define N        2         /*liczba planet*/
#define LP       11        /*liczba parametrow*/
#define M        38        /*liczba pomiarow*/ 
#define B        0.0001    /*dokladnosc rozwiazania r-nia Keplera*/
#define KROK     0.001     /*podzial siatki argumentow*/
#define ZAKRES   0.05      /*zakres zmian parametrow*/
#define PI       3.141592653589793115997963468544185161590576172 /*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);


main()
{
 int i,j;
 float Fn;
 float V0,K[N],P[N],e[N],w[N],tau[N];
 float s[LP],s0[LP],t[M],v[M],d[M];
 FILE *plik, *plik1, *plik2, *plik3, *plik4, *plik5, *plik6, *plik7, *plik8;
 FILE *plik9, *plik10, *plik11, *wynik;
 char *nazwa;
 float param[LP]={1,1,DZIEN,1,PI/180,DZIEN,1,DZIEN,1,PI/180,DZIEN};



/*wczytanie danych*/

 nazwa="HD114729.dat";
 plik=fopen(nazwa,"r");
 
 for (i=0; i<=M-1; ++i)
 {
  fscanf(plik,"%f%f%f",&t[i],&v[i],&d[i]);
  t[i]=t[i]*DZIEN;
 }
 
 fclose(plik);
 
 
/*parametry poczatkowe*/

 V0=-5.32061;
 K[0]=15.29573;
 P[0]=1088.65051*DZIEN;
 e[0]=0.205788;
 w[0]=75.49266*PI/180;
 tau[0]=509.81573*DZIEN;
 K[1]=5.69828;
 P[1]=10.52720*DZIEN;
 e[1]=0.016552;
 w[1]=38.30322*PI/180;
 tau[1]=100.91705*DZIEN;

 nazwa="HD114729_chi2_V0_1.dat";
 plik1=fopen(nazwa,"w");
 nazwa="HD114729_chi2_K1_1.dat";
 plik2=fopen(nazwa,"w");
 nazwa="HD114729_chi2_P1_1.dat";
 plik3=fopen(nazwa,"w");
 nazwa="HD114729_chi2_e1_1.dat";
 plik4=fopen(nazwa,"w");
 nazwa="HD114729_chi2_w1_1.dat";
 plik5=fopen(nazwa,"w");
 nazwa="HD114729_chi2_tau1_1.dat";
 plik6=fopen(nazwa,"w");
 nazwa="HD114729_chi2_K2_1.dat";
 plik7=fopen(nazwa,"w");
 nazwa="HD114729_chi2_P2_1.dat";
 plik8=fopen(nazwa,"w");
 nazwa="HD114729_chi2_e2_1.dat";
 plik9=fopen(nazwa,"w");
 nazwa="HD114729_chi2_w2_1.dat";
 plik10=fopen(nazwa,"w");
 nazwa="HD114729_chi2_tau2_1.dat";
 plik11=fopen(nazwa,"w");
 

 /* glowna czesc programu*/
 
 s0[0]=V0;
 for (j=0; j<=N-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 (i=0; i<=LP-1; ++i)
{
 for (j=0; j<=LP-1; ++j)
 {
  s[j]=s0[j];
 }
 s[i]=(1-ZAKRES)*s0[i];
 
 if (i==0) wynik=plik1;
  
 if (i==1) wynik=plik2;
  
 if (i==2) wynik=plik3;
  
 if (i==3) wynik=plik4;
  
 if (i==4) wynik=plik5;
  
 if (i==5) wynik=plik6;
  
 if (i==6) wynik=plik7;
  
 if (i==7) wynik=plik8;
  
 if (i==8) wynik=plik9;
  
 if (i==9) wynik=plik10;
  
 if (i==10) wynik=plik11;
 
 if (s0[i]<0)
 {
  while (s[i]>((1+ZAKRES)*s0[i]))
  {
   Fn=f(s,t,v,d,M,LP,B,N,PI);
   fprintf(wynik,"%10.5f%10.5f\n",s[i]/param[i],Fn);
   s[i]=s[i]+KROK*s0[i];
  }
 }
 else
 {
  while (s[i]<((1+ZAKRES)*s0[i]))
  {
   Fn=f(s,t,v,d,M,LP,B,N,PI);
   fprintf(wynik,"%10.5f%10.5f\n",s[i]/param[i],Fn);
   s[i]=s[i]+KROK*s0[i];
  }
 }
}

fclose(plik1);
fclose(plik2);
fclose(plik3);
fclose(plik4);
fclose(plik5);
fclose(plik6);
fclose(plik7);
fclose(plik8);
fclose(plik9);
fclose(plik10);
fclose(plik11);

return 0;

}

/*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;
} 


-- CezaryMigaszewski - 26 Feb 2005

This topic: Main > TWikiUsers > CezaryMigaszewski > FunkcjaDopasowania
Topic revision: 05 Mar 2005, CezaryMigaszewski
 
This site is powered by FoswikiCopyright © CC-BY-SA by the contributing authors. All material on this collaboration platform is copyrighted under CC-BY-SA by the contributing authors unless otherwise noted.
Ideas, requests, problems regarding Foswiki? Send feedback