```\ Program solution to PHYS 551 Term Project 2000

FALSE [IF]

Find the initial angles for a cannon ball fired from
a cannon located at 43 degrees N. latitude, with
elevation theta and angle (relative to true North) phi,
to reach a target 2000 m due East and 100 m higher,
within 0.5 m. Also determine the sensitivity of the
result to variations of the muzzle velocity and the
aiming angles.

My solution:

1. The 6 equations of motion are

dv_x/dt + D*v*v_x - 2*w*(S*v_y - C*v_z) = 0

dv_y/dt + D*v*v_y + 2*w*v_x*S = 0

dv_z/dt + D*v*v_z + g - 2*w*v_x*C = 0

dx/dt = v_x, etc.

where

D = .5*c*A*rho(z)/M ,

v = sqrt(v_x^2 + v_y^2 + v_z^2) ,

rho(z) = rho_0 * exp(-z/L)

is the atmospheric density varying with altitude z ,

S = sin(chi), C = cos(chi)
(chi is the conventional latitude, measured from
the Equator)

and w is the Earth's angular velocity.

2. Change independent variable from time to x-coordinate.

3. Compute the constants and initial conditions:

v = 180 m/s

v_x = v * sin(theta) * cos(phi)
v_y = v * sin(theta) * sin(phi)
v_z = v * cos(theta)

(x,y,z) = (0,0,0)

S = sin(43 deg)     C = cos(43 deg)

4. Solve by 2nd-order Runge-Kutta, with stepsize <= 0.25 m

5. Find minimum of |vec(x) - vec(x_final)|^2
(the method used below is quadratic minimization,
applied recursively, first to elevation theta and
then to aiming angle phi).

Parameters: M = 20 kg           \ projectile mass
c = 0.3             \ shape factor
rho_Pb = 11.35      \ density of Pb
A = 0.0176 m^2      \ cross-section
g = 9.80665 m/s^2   \ accel of gravity
x_fin = 2000 m      \ range
y_fin = 0           \ N-S offset
zfin = 100 m        \ target height

[THEN]

MARKER  -proj

INCLUDE ftran111.f
INCLUDE flocals.f

\ INCLUDE ansfalsi.f
\ INCLUDE ctran3.f

\$" fvalue"  FIND  NIP  0=  [IF]  INCLUDE fvalue.f  [THEN]

: fvariables     0 DO  FVARIABLE  LOOP  ;

7 fvariables x xstep y dy z dz t            \ coordinates
6 fvariables Vx Vy Vz kx ky kz              \ velocities
7 fvariables V0 V Vp dV Vxp Vyp  Vzp        \ more velocities
7 fvariables D1 zscale g M rho Area Shape   \ parameters
5 fvariables drag  half  wx2  invVx invVxp  \ misc. variables
2 fvariables wx2S  wx2C                     \ latitude params
3 fvariables xfin  yfin  zfin               \ impact point

1e0 F2/  half F!                            \ multiplier by 0.5

\$" FPI"   FIND NIP  0= [IF]  PI  FCONSTANT  FPI  [THEN]
: radians   ( f: degrees -- radians)  FPI  F*  180e0  F/  ;

: set_params
f" M = 20"          f" Shape = 0.3"
f" Area = 0.0176"   f" rho = 0.029/22.4E-3"
f" g = 9.80665"

f" wx2 = 4 * 3.1415927 / (24 * 3600) "
f" wx2S = wx2 * sin(radians(43)) "
f" wx2C = wx2 * cos(radians(43)) "

f" zscale = 8.31662 * 293 / (0.029 * g)"

f" D1 = 0.5 * Shape * Area * rho / M"

f" xfin = 2000"     f" yfin = 0"    f" zfin = 100"

f" V0 = 180"
;

: D     ( f: z -- D)    \ f" D1 * exp(-z/zscale) "
zscale F@ F/       ( f: -- u)
( FNEGATE   FEXP   D1 F@   F* )
FDUP   half F@ F*      ( f: -- u u/2)
1e0  F-  F*  1e0  F+   ( f: -- 1-u*[1-u/2])
D1  F@  F*
\ approximation good to 2e-3 at z = 2 Km
;

: first_step   \ 1st R-K step
f" drag = D(z) * V"
f" invVx = 1/Vx"
f" kx = -xstep * (drag - invVx * (wx2S*Vy - wx2C*Vz)) "
f" ky = -xstep * (drag * Vy * invVx + wx2S)"
f" kz = -xstep * ( (drag * Vz + g)*invVx - wx2C )"

f" Vxp = Vx + kx"
f" Vyp = Vy + ky"
f" Vzp = Vz + kz"

f" Vp = V + (Vx*kx + Vy*ky + Vz*kz) / V"
;

: second_step   \ 2nd R-K step
f" dy = xstep * Vy * invVx"
f" dz = xstep * Vz * invVx"

f" drag = D(z+dz) * Vp"

f" invVxp = 1/Vxp"

f" kx = -xstep * (drag - invVxp * (wx2S*Vyp - wx2C*Vzp) ) "
f" ky = -xstep * (drag * Vyp * invVxp + wx2S)"
f" kz = -xstep * ( (drag * Vzp + g)*invVxp - wx2C )"

\ update velocities
f" Vx = half * ( Vxp + Vx + kx )"
f" Vy = half * ( Vyp + Vy + ky )"
f" Vz = half * ( Vzp + Vz + kz )"
f" V  = half * ( V + Vp + (Vxp*kx + Vyp*ky + Vzp*kz) / Vp )"

\ update coords
f" z = z + half * ( dz + xstep * Vz/Vx )"
f" y = y + half * ( dy + xstep * Vy/Vx )"
f" t = t + xstep/Vx"
f" x = x + xstep"
;

: new_point
first_step   second_step
;

set_params

: initialize    ( f: theta phi -- )     \ initialize variables
frame| aa bb |                      \ aa = phi, bb = theta
\ theta is elevation
f" V = V0"

f" Vz = V * sin(radians(bb)) "

f" x=0"
f" y=0"
f" z=0"
f" t=0"
|frame
;

: fire     ( f: theta phi -- )  \ integrate diffeqs until impact
initialize
new_point
BEGIN  z F@  zfin F@  F>    Vz  F@ F0>  OR
WHILE  new_point  REPEAT
;

: f1    ( f: theta phi -- residual)
fire                        \ shoot cannon
f" (x-xfin)^2 + (y-yfin)^2 + (z-zfin)^2 "
;

: .vec   CR   x F@ F.   y F@ F.   z F@ F.  ;

: .vel   CR  Vx F@ F.  Vy F@ F.  Vz F@ F.  ;

2 fvariables  theta phi

: f2    ( f: theta --)   phi  F@   f1  ;

: f3    ( f: phi --)  theta F@   FSWAP   f1  ;

: mintheta  ( f: theta_min theta_max phi --)  \ minimize w.r.t. theta
frame| aa bb cc |    \ aa = phi, bb = theta_max...
f" phi = aa"
use( f2 cc F@  bb F@  )quadmin
|frame
;

: minphi    ( f: phi_min  phi_max theta --)   \ minimize w.r.t. phi
frame| aa bb cc |
f" theta = aa "
use( f3 cc F@  bb F@  )quadmin
|frame
;
```