20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
|
# File 'ext/ballistics/ext.c', line 20
VALUE method_trajectory(VALUE self,
VALUE drag_function,
VALUE drag_coefficient,
VALUE velocity,
VALUE sight_height,
VALUE shooting_angle,
VALUE zero_angle,
VALUE wind_speed,
VALUE wind_angle,
VALUE max_range,
VALUE interval) {
/* cast ruby variables */
int DragFunction = FIX2INT(drag_function);
double DragCoefficient = NUM2DBL(drag_coefficient);
double Vi = NUM2DBL(velocity);
double SightHeight = NUM2DBL(sight_height);
double ShootingAngle = NUM2DBL(shooting_angle);
double ZAngle = NUM2DBL(zero_angle);
double WindSpeed = NUM2DBL(wind_speed);
double WindAngle = NUM2DBL(wind_angle);
int MaxRange = FIX2INT(max_range);
int Interval = FIX2INT(interval);
/* create ruby objects */
VALUE result_array = rb_ary_new2(MaxRange);
double t=0;
double dt=0.5/Vi;
double v=0;
double vx=0, vx1=0, vy=0, vy1=0;
double dv=0, dvx=0, dvy=0;
double x=0, y=0;
double headwind=HeadWind(WindSpeed, WindAngle);
double crosswind=CrossWind(WindSpeed, WindAngle);
double Gy=GRAVITY*cos(DegtoRad((ShootingAngle + ZAngle)));
double Gx=GRAVITY*sin(DegtoRad((ShootingAngle + ZAngle)));
vx=Vi*cos(DegtoRad(ZAngle));
vy=Vi*sin(DegtoRad(ZAngle));
y=-SightHeight/12;
int n=0;
for (t=0;;t=t+dt){
vx1=vx, vy1=vy;
v=pow(pow(vx,2)+pow(vy,2),0.5);
dt=0.5/v;
// Compute acceleration using the drag function retardation
dv = retard(DragFunction,DragCoefficient,v+headwind);
dvx = -(vx/v)*dv;
dvy = -(vy/v)*dv;
// Compute velocity, including the resolved gravity vectors.
vx=vx + dt*dvx + dt*Gx;
vy=vy + dt*dvy + dt*Gy;
int yards = (x/3);
if (yards>=n){
if (yards % Interval == 0){
VALUE entry = rb_hash_new();
double windage_value = Windage(crosswind,Vi,x,t+dt);
double moa_windage_value = windage_value / ((yards / 100.0) * 1.0465);
rb_hash_aset(entry,
rb_str_new2("range"),
rb_float_new((int)(yards)));
rb_hash_aset(entry,
rb_str_new2("path"),
rb_float_new(y*12));
rb_hash_aset(entry,
rb_str_new2("moa_correction"),
rb_float_new(-RadtoMOA(atan(y/x))));
rb_hash_aset(entry,
rb_str_new2("time"),
rb_float_new(t+dt));
rb_hash_aset(entry,
rb_str_new2("windage"),
rb_float_new(windage_value));
rb_hash_aset(entry,
rb_str_new2("moa_windage"),
rb_float_new(moa_windage_value));
rb_hash_aset(entry,
rb_str_new2("velocity"),
rb_float_new(v));
rb_ary_push(result_array, entry);
}
n++;
}
// Compute position based on average velocity.
x=x+dt*(vx+vx1)/2;
y=y+dt*(vy+vy1)/2;
if (fabs(vy)>fabs(3*vx)) break;
if (n>=MaxRange+1) break;
}
return result_array;
}
|