Program sluzy do rysowania krzywej modelowej predkosci radialnej gwiazdy, wywolanej istnieniem planet wokol tej gwiazdy. Wszystkie parametry trzeba wpisac w kodzie ( program o nic nas nie pyta po uruchomieniu ). Nalezy ustawic liczbe planet N ( w przykladzie ponizej sa dwie planety ), poczatkowy czas t, liczbe krokow czasowych LK, dlugosc kroku h. Nalezy takze podac parametry orbitalne V0,K,P,e,w,tau(w kodzie podalem przykladowe parametry i wyjanilem co oznaczaja). Nalezy takze podac nazwe piku wyjsciowego, do ktorego zostana zapisane punkty krzywej predkosci radialnej. W programie mozna zmieniac czas poczatkowy krzywej, liczbe i dlugosc kroku czasowego.
/ *modelowa_krzywa.c* /
#include<stdio.h>
#include<math.h>
#define DZIEN 86400 / *liczba sekund w 1 dniu* /
#define N 2 / *liczba planet* /
#define B 0.0001 / *dokladnosc rozwiazania r-nia Keplera* /
#define LK 25000 / *liczba krokow* /
#define PI 3.1415926536 / *liczba pi* /
float NR(float P, float tau, float t, float e, float b, float pi);
main()
{
int i,k;
float h,t,Vr,niu,E;
float V0,K[N],P[N],e[N],w[N],tau[N];
FILE *wynik;
char *nazwa;
t=400*DZIEN; /* poczatkowy czas */
h=0.1*DZIEN; /* dlugosc kroku */
/*parametry orbitalne - tu dokonujemy zmian */
V0=-5.5; /* predkosc srodka masy */
K[0]=15; /* ampituda predkosci radialnej wywolana 1. planeta */
P[0]=1080*DZIEN; /* okres 1. planety */
e[0]=0.19; /* ekscentycznosc orbity 1. planety */
w[0]=81*PI/180; /* dlugosc perycentrum orbity 1. planety */
tau[0]=539*DZIEN; /* czas przejscia przez pericentrum 1. plnety */
K[1]=6; /* podobnie dla kolejnych planet */
P[1]=10.5262*DZIEN; /* -- */
e[1]=0.016; /* -- */
w[1]=53*PI/180; /* -- */
tau[1]=101.4*DZIEN; /* -- */
nazwa="HD114729_modelowa.dat"; /* nazwa pliku wyjsciowego */
wynik=fopen(nazwa,"w");
for (k=1; k<=LK; ++k)
{
Vr=V0;
for (i=0; i<=N-1; ++i)
{
E=NR(P[i],tau[i],t,e[i],B,PI);
niu=2*atan(sqrt((1+e[i])/(1-e[i]))*tan(E/2));
Vr=Vr+K[i]*(cos(w[i]+niu)+e[i]*cos(w[i]));
}
fprintf(wynik,"%10.4f%10.2f%10.2f\n",t/DZIEN,Vr,1);
t=t+h;
}
fclose(wynik);
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;
}
--
CezaryMigaszewski - 26 Feb 2005
wykresy
Commands to make figures - tac is cat backwards, rób
info tac lub
info tail lub
info awk dla więcej info. Też
graph jest części
plotutils ComputerLanguages#plotutils_GPL_plotting_library, możesz używać albo na linii komenda, albo w programie jak biblioteki. Też możesz zrobić
graph -T ps dla plik postscript, lub
graph -T X dla okno na ekranu, itd...
- tac HD114729_modelowa.dat |tail -500 | awk '{print $1 " " $2}' | graph -T png > HD114729_500.png
- HD114729_500.png:
- tac HD114729_modelowa.dat |tail -5000 | awk '{print $1 " " $2}' | graph -T png > HD114729_5000.png
- HD114729_5000.png:
- cat HD114729_modelowa.dat | awk '{print $1 " " $2}' | graph -T png > HD114729_all.png
- HD114729_all.png:
syntaks
Z zwykly kompilator (gcc-2.95.4), syntaks
/ * komentarz tu * / z spacjami nie jest akceptowany i daje dziwny błędy w kompilacji. Standardowy syntaks
/* komentarz tu */ kompiluje fajnie. Chyba masz kompilator, który akceptuje pierwsza wersję.