Skip to content

Commit 3afc04e

Browse files
author
Cosmicstring
committed
Uploading all the code
1 parent c8fdc6e commit 3afc04e

20 files changed

+2323
-1
lines changed

‎.gitignore

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
*.txt
2+
*.in
3+
*.o
4+
*.out
5+
*.aei

‎Binning_Meteoroids.cpp

Lines changed: 161 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,161 @@
1+
#include <stdio.h>
2+
#include <stdlib.h>
3+
#include <math.h>
4+
#include <string.h>
5+
6+
#define MaxN 2001
7+
#define resolution 1024
8+
#define DEBUG 1
9+
10+
11+
long A[MaxN*MaxN];
12+
long P1[MaxN*MaxN], P2[MaxN*MaxN];
13+
long Test[MaxN];
14+
15+
FILE* SaSOHO;
16+
FILE* Presek1;
17+
FILE* Presek2;
18+
FILE* TestingOutput;
19+
FILE* DebugText;
20+
FILE* Presek1Ispis;
21+
FILE* Presek2Ispis;
22+
23+
int main(){
24+
25+
Presek1 = fopen("Trenutak1.txt", "r");
26+
Presek2 = fopen("Trenutak2.txt", "r");
27+
Presek1Ispis = fopen("Presek1Ispis.txt", "w");
28+
Presek2Ispis = fopen("Presek2Ispis.txt", "w");
29+
SaSOHO = fopen("SaSOHO.txt", "r");
30+
TestingOutput = fopen("TestingOutput_Beta06.txt", "w"); /**Od najvece dimenzije do najmanje 3-0 {0.05, 0.5, 0.6, 1}*/
31+
DebugText = fopen("Debug.txt", "w");
32+
double foo, xi, yi, Testx, Testy, x,y, MinX, MinY, MaxX, MaxY, binstep, binstepP1;
33+
long i,j;
34+
char input[256];
35+
36+
MaxX = 17; MaxY = 17;
37+
MinX = -17; MinY = -17; /** [stepeni] Vrednost celog vidnog polja sa Suncem u Centru, Sekanin rad za velicinu piksela(Uzeo sam kameru C3)*/
38+
memset(Test, sizeof(Test), 0);
39+
memset(A, sizeof(A), 0);
40+
binstep = (MaxX-MinX)/resolution; /**Kamera C3 ima 1024*1024 piksela*/
41+
42+
if(DEBUG) printf("%lf\n", binstep);
43+
44+
/*for(i=0; i<10; i++){
45+
/*fgets(input, sizeof input, Testing);
46+
sscanf(input, "%lf %lf", &Testx, &Testy);
47+
scanf("%lf%lf", &Testx, &Testy);
48+
printf("%lf\t%ld\t%lf\t%ld\n", fabs(Testx)/binstep, lround(fabs(Testx)/binstep), fabs(Testy)/binstep, lround(fabs(Testy)/binstep));
49+
printf("1: %ld\n", lround(fabs(Testx)/binstep) + lround(fabs(Testy)/binstep)*10);
50+
Test[lround(fabs(Testx)/binstep) + lround(fabs(Testy)/binstep)*10]++;
51+
}
52+
53+
for(i=0; i<10; i++)
54+
for(j=0; j<10; j++){
55+
if(j==9) printf("%ld\n", Test[i+j*10]);
56+
else printf("%ld ", Test[i+j*10]);
57+
}*/
58+
59+
while(!feof(SaSOHO)){
60+
fgets(input, sizeof input, SaSOHO);
61+
sscanf(input, "%lf %lf", &x, &y);
62+
// if(DEBUG) printf("Ok: %ld %ld %ld\n", A[lround(fabs(x)/binstep) + lround((fabs(y)/binstep))*2000], lround(fabs(x)/binstep) + lround((fabs(y)/binstep))*2000, i);
63+
if(x<0) xi = MaxX - fabs(x);
64+
else
65+
if(x>0) xi = MaxX+x;
66+
else xi = x;
67+
if(y<0) yi = MaxY+fabs(y);
68+
else
69+
if(y>0) yi = MaxY - y;
70+
else yi = y;
71+
A[lround(xi/binstep) + lround(yi/binstep)*resolution]++;
72+
// if(DEBUG) printf("Ok: %ld %ld %ld\n", A[lround(fabs(x)/binstep) + lround((fabs(y)/binstep))*2000], lround(fabs(x)/binstep) + lround((fabs(y)/binstep))*2000, i);
73+
}
74+
75+
76+
for(i=0; i<resolution*resolution; i++){
77+
if(i%resolution==0) fprintf(TestingOutput, "%ld\n", A[i]);
78+
else fprintf(TestingOutput, "%ld\t", A[i]);
79+
}
80+
81+
/** Presek 1 */
82+
83+
memset(input, sizeof input, 0);
84+
memset(P1, sizeof P1, 0);
85+
double MaxXP1, MinXP1, MaxYP1, MinYP1;
86+
double binstepP1X, binstepP1Y;
87+
88+
MaxXP1 = -100; MaxYP1 = -100;
89+
MinXP1 = 100; MinYP1 = 100;
90+
91+
while(!feof(Presek1)){
92+
fgets(input, sizeof input, Presek1);
93+
sscanf(input, "%lf %lf %lf", &foo, &x, &y);
94+
if(MaxXP1<fabs(x)) MaxXP1 = fabs(x);
95+
if(MaxYP1<fabs(y)) MaxYP1 = fabs(y);
96+
if(MinXP1>fabs(x)) MinXP1 = fabs(x);
97+
if(MinYP1>fabs(y)) MinYP1 = fabs(y);
98+
}
99+
100+
binstepP1X = (fabs(MaxXP1)-fabs(MinXP1))/resolution;
101+
binstepP1Y = (fabs(MaxYP1)-fabs(MinYP1))/resolution;
102+
if(DEBUG) printf("%.10lf %lf %lf %lf %lf %lf %lf %lf\n", binstepP1X, fabs(MaxXP1)-fabs(MinXP1), (fabs(MaxYP1) - fabs(MinYP1))/binstepP1X, (fabs(MaxXP1)-fabs(MinXP1))/binstepP1Y, MaxXP1, MaxYP1, MinXP1, MinYP1);
103+
rewind(Presek1);
104+
105+
106+
i=0;
107+
while(!feof(Presek1)){
108+
fgets(input, sizeof input, Presek1);
109+
sscanf(input, "%lf %lf %lf", &foo, &x, &y);
110+
P1[lround((fabs(x)-fabs(MinXP1))/binstepP1X) + lround((fabs(y)-fabs(MinYP1))/binstepP1Y)*resolution]++;
111+
}
112+
113+
for(i=0; i<resolution*resolution; i++)
114+
if(i%resolution==0) fprintf(Presek1Ispis, "%ld\n", P1[i]);
115+
else fprintf(Presek1Ispis, "%ld\t", P1[i]);
116+
117+
/** Presek 2*/
118+
119+
memset(input, sizeof input, 0);
120+
memset(P2, sizeof P2, 0);
121+
double MaxXP2, MinXP2, MaxYP2, MinYP2;
122+
double binstepP2X, binstepP2Y;
123+
124+
MaxXP2 = -100; MaxYP2 = -100;
125+
MinXP2 = 100; MinYP2 = 100;
126+
127+
while(!feof(Presek2)){
128+
fgets(input, sizeof input, Presek2);
129+
sscanf(input, "%lf %lf %lf", &foo, &x, &y);
130+
if(MaxXP2<fabs(x)) MaxXP2 = fabs(x);
131+
if(MaxYP2<fabs(y)) MaxYP2 = fabs(y);
132+
if(MinXP2>fabs(x)) MinXP2 = fabs(x);
133+
if(MinYP2>fabs(y)) MinYP2 = fabs(y);
134+
}
135+
136+
if(DEBUG) printf("%lf %lf %lf %lf\n", MaxXP2, MinXP2, MaxYP2, MinYP2);
137+
138+
binstepP2X = (fabs(MaxXP2)-fabs(MinXP2))/resolution;
139+
binstepP2Y = (fabs(MaxYP2)-fabs(MinYP2))/resolution;
140+
if(DEBUG) printf("%.10lf %lf %lf %lf %lf %lf %lf %lf\n", binstepP2X, fabs(MaxXP2)-fabs(MinXP2), (fabs(MaxYP2) - fabs(MinYP2))/binstepP2X, (fabs(MaxXP2)-fabs(MinXP2))/binstepP2Y, MaxXP2, MaxYP2, MinXP2, MinYP2);
141+
142+
rewind(Presek2);
143+
144+
while(!feof(Presek2)){
145+
fgets(input, sizeof input, Presek2);
146+
sscanf(input, "%lf %lf %lf", &foo, &x, &y);
147+
P2[lround((fabs(x)-fabs(MinXP2))/binstepP2X) + lround((fabs(y)-fabs(MinYP2))/binstepP2Y)*resolution]++;
148+
}
149+
150+
for(i=0; i<resolution*resolution; i++)
151+
if(i%resolution==0) fprintf(Presek2Ispis, "%ld\n", P2[i]);
152+
else fprintf(Presek2Ispis, "%ld\t", P2[i]);
153+
154+
fclose(SaSOHO);
155+
fclose(Presek1Ispis);
156+
fclose(Presek1);
157+
fclose(Presek2);
158+
fclose(DebugText);
159+
fclose(TestingOutput);
160+
return 0;
161+
}

‎Binning_SOHO.cpp

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
#include <iostream>
2+
#include <cstdio>
3+
#include <cstdlib>
4+
#include <cmath>
5+
#include <cstring>
6+
7+
/**
8+
9+
THIS IS A CODE WRITTEN FOR MAKING A BINNED IMAGE FROM THE XYZ DATA OF DUST PARTICLES
10+
EJECTED FROM THE COMET SURFACE. THE IMAGE IS LATER ON FORWARDED TO MATLAB CODES WHICH
11+
PROCESS THEM FURTHER.
12+
13+
**/
14+
15+
16+
#define MaxN 2001
17+
#define DEBUG 0
18+
19+
20+
long A[MaxN*MaxN];
21+
long Test[MaxN];
22+
23+
FILE* SaSOHO;
24+
FILE* Testing;
25+
FILE* TestingOutput;
26+
FILE* DebugText;
27+
28+
int main(){
29+
30+
SaSOHO = fopen("SaSOHO.txt", "r");
31+
Testing = fopen("Test.txt", "r");
32+
TestingOutput = fopen("TestingOutput.txt", "w");
33+
DebugText = fopen("Debug.txt", "w");
34+
double Testx, Testy, x,y, MinX, MinY, MaxX, MaxY, binstep;
35+
long i,j;
36+
char input[256];
37+
38+
MaxX = 0.2; MaxY = 0.2;
39+
MinX = -0.2; MinY = -0.2; /** Vrednost celog vidnog polja sa Suncem u Centru */
40+
memset(Test, sizeof(Test), 0);
41+
memset(A, sizeof(A), 0);
42+
binstep = (MaxX-MinX)/2000;
43+
44+
if(DEBUG) printf("%lf\n", binstep);
45+
46+
/*
47+
TESTING ROUTINE
48+
49+
for(i=0; i<10; i++){
50+
/*fgets(input, sizeof input, Testing);
51+
sscanf(input, "%lf %lf", &Testx, &Testy);
52+
scanf("%lf%lf", &Testx, &Testy);
53+
printf("%lf\t%ld\t%lf\t%ld\n", fabs(Testx)/binstep, lround(fabs(Testx)/binstep), fabs(Testy)/binstep, lround(fabs(Testy)/binstep));
54+
printf("1: %ld\n", lround(fabs(Testx)/binstep) + lround(fabs(Testy)/binstep)*10);
55+
Test[lround(fabs(Testx)/binstep) + lround(fabs(Testy)/binstep)*10]++;
56+
}
57+
58+
for(i=0; i<10; i++)
59+
for(j=0; j<10; j++){
60+
if(j==9) printf("%ld\n", Test[i+j*10]);
61+
else printf("%ld ", Test[i+j*10]);
62+
}*/
63+
i=0;
64+
while(!feof(SaSOHO)){
65+
fgets(input, sizeof input, SaSOHO);
66+
sscanf(input, "%lf %lf", &x, &y);
67+
if(DEBUG) printf("Ok: %ld %ld %ld\n", A[lround(fabs(x)/binstep) + lround((fabs(y)/binstep))*2000], lround(fabs(x)/binstep) + lround((fabs(y)/binstep))*2000, i);
68+
A[lround(fabs(x)/binstep) + lround((fabs(y)/binstep))*2000]++;
69+
if(DEBUG) printf("Ok: %ld %ld %ld\n", A[lround(fabs(x)/binstep) + lround((fabs(y)/binstep))*2000], lround(fabs(x)/binstep) + lround((fabs(y)/binstep))*2000, i);
70+
i++;
71+
}
72+
73+
for(i=0; i<2000*2000; i++){
74+
if((i%2000==0)) fprintf(TestingOutput, "%ld\n", A[i]);
75+
else fprintf(TestingOutput, "%ld\t", A[i]);
76+
}
77+
78+
if(DEBUG) printf("OkIzaso %ld\n", A[1206304]);
79+
/*
80+
for(i=0; i<2000*2000; i++){
81+
if(i%2000==0) fprintf(TestingOutput, "%ld\n", A[i]);
82+
else fprintf(TestingOutput, "%ld\t", A[i]);
83+
}*/
84+
85+
if(DEBUG) printf("Nema SegFault\n");
86+
87+
fclose(SaSOHO);
88+
fclose(Testing);
89+
fclose(DebugText);
90+
fclose(TestingOutput);
91+
return 0;
92+
}

‎Comet_Particle_Ejection.cpp

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
#include "Comet_Particle_Ejection.h"
2+
3+
double Z1,ZBB,TBB,PvBB,PvZ1,T1;
4+
double molmasavode = (18.02*1.661*1e-27);
5+
double BolzConst = 1.3806488e-23;
6+
// Just to eject one particle at a time which is according to Whipple formulae
7+
int ni=0;
8+
9+
//To calculate the total production, i.e. total ejected mass
10+
double Zuu = 0;
11+
12+
int main()
13+
{
14+
//Declaration of I/O files
15+
FILE* Provera = fopen("Provera.txt","w");
16+
FILE* CheckInput = fopen("CheckInput.txt","w");
17+
FILE* Vecheck = fopen("Vecheck.txt","w");
18+
FILE* BrzinaiPolozaj_MercuryCestice = fopen("beta0_5_small.in","w");
19+
FILE* Churyumov_Gerasimenko = fopen("67PComet.aei", "r");
20+
FILE* Comet_Cartesian = fopen("Cartesian_Comet.txt", "w");
21+
22+
OrbitalElements p, o1,o;
23+
PhaseState p1,p2,pcomet;
24+
25+
comet Comet_CG;
26+
27+
vector r,v;
28+
29+
int i,k,kt=1,j;
30+
int MeteorNum = 0;
31+
32+
char input[200];
33+
34+
double Te = 365.256;
35+
double ae = 149.6e9;
36+
double K= Te/sqrt(ae*ae*ae);
37+
double N= 5120;
38+
double u1,u,u0,eps,mi,m;
39+
double f0,f1,f2,f3,d1,d2,d3;
40+
double TA,T,TA0,ta,ta0;
41+
double t, t0;
42+
double rp,x3,y3,z3,x2,y2,z2,x1,y1,z1,x,y,z;
43+
double lambda,beta;
44+
double Q;
45+
double timeobs = 365.256;
46+
double h = 0.01;
47+
double dt1, dt2;
48+
double Ugao,Ri;
49+
50+
double Time;
51+
52+
srand(time(NULL));
53+
54+
/*From paper: A homogeneous nucleus for comet 67P/Churyumov–Gerasimenko from its gravity field*/
55+
56+
Comet_CG.R = 2.65;
57+
Comet_CG.ro = 0.553e12;
58+
Comet_CG.Mass = 9.982e9;
59+
60+
/*From JPL*/
61+
p.e = 0.64093009801435;
62+
p.a = 3.462505635805988; /**[AU]*/
63+
p.i = ConvertToRadians(7.040168070284028);
64+
p.argperi = ConvertToRadians(12.80065127074631);
65+
p.longnode = ConvertToRadians(50.13500663540592);
66+
67+
cartesian(2.959122082855911e-4, p, &pcomet);
68+
69+
fprintf(Comet_Cartesian, "%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\t%.10lf\n", pcomet.x, pcomet.y, pcomet.z, pcomet.xd, pcomet.yd, pcomet.zd);
70+
71+
/*Skip first 4 lines*/
72+
for(int i=1; i<5; i++) fgets(input, sizeof input, Churyumov_Gerasimenko);
73+
74+
//Set up the file for mercury
75+
fprintf(BrzinaiPolozaj_MercuryCestice, ")O+_06 Small-body initial data (WARNING: Do not delete this line!!)\nstyle (Cartesian, Asteroidal, Cometary) = Ast\n)---------------------------------------------------------------------\n");
76+
77+
while(!feof(Churyumov_Gerasimenko)){
78+
79+
fgets(input, sizeof input, Churyumov_Gerasimenko);
80+
if(DEBUG) printf("Ok1\n");
81+
82+
sscanf(input, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", &Time, &r.x, &r.y, &r.z, &v.x, &v.y, &v.z, &Comet_CG.a, &Comet_CG.e, &Comet_CG.i, &Comet_CG.argperi, &Comet_CG.longnode, &Comet_CG.meananom);
83+
if(DEBUG) printf("Ok2\n");
84+
85+
fprintf(CheckInput, "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", Time, r.x, r.y, r.z, v.x, v.y, v.z, Comet_CG.a, Comet_CG.e, Comet_CG.i, Comet_CG.argperi, Comet_CG.longnode, Comet_CG.meananom);
86+
87+
/*Vracam u km i km/d posto se sa tim radi u GiveMeParticles*/
88+
double AUkm = AU * 1e-3;
89+
r.x = r.x * AUkm;
90+
r.y = r.y * AUkm;
91+
r.z = r.z * AUkm;
92+
93+
v.x = v.x * AUkm;
94+
v.y = v.y * AUkm;
95+
v.z = v.z * AUkm;
96+
97+
//Time in days so it can be added to Epoch of the comet
98+
GiveMeParticles(BrzinaiPolozaj_MercuryCestice, Vecheck, Provera, MeteorNum, Zuu, ni, Comet_CG.Mass, Comet_CG.R, Time*365.25, p, r, v);
99+
100+
}
101+
102+
if(DEBUG) printf("%.10lf\n", Comet_CG.Mass);
103+
104+
fclose(Provera);
105+
fclose(Vecheck);
106+
fclose(CheckInput);
107+
fclose(Churyumov_Gerasimenko);
108+
fclose(BrzinaiPolozaj_MercuryCestice);
109+
fclose(Comet_Cartesian);
110+
111+
return 0;
112+
}

0 commit comments

Comments
 (0)