\ 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 quadmin1.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" Vx = V * cos(radians(aa)) * cos(radians(bb)) "
    f" Vy = V*sin(radians(aa)) * cos(radians(bb)) "
    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
;