\ 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 ;