You are here: Foswiki>Main Web>TWikiUsers>CezaryMigaszewski>FunkcjaDopasowania (revision 1)EditAttach
Omowienie programu wkrotce



/*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
Edit | Attach | Print version | History: r2 < r1 | Backlinks | View wiki text | Edit WikiText | More topic actions...
Topic revision: r1 - 26 Feb 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