General remarks first, applicable to all assignments. Please include your thoroughly documented/commented code and description of methods and results in a Jupyter notebook submitted via Quercus to TA Fergus H. Every important line needs to be commented to show that you correctly understand the algorithm (do not comment "z gets multiplied by y". We can read that. We what to know why you do certain things. DEADLINE prolonged 2 days to: 28 November 2019 in the evening WARNING: We do run a software that is very good at spotting similarities of your homework to somebody else's code. Discussing methods with others is ok, but implementing (writing actual code) with someone or copying code's details is plagiarism. Problem 1 [5 pt] - Integrate diffraction pattern _________________________________________________ A laser beam of coherent red light with wavelength 0.6 um (micrometers) passes at right angle though a single slit 30 um wide. Compute the 1-dimensional profile of intensity of light on a screen placed z = 1 m away from the slit, covering distance y=-6 to y=+6 cm off beam axis. Normalize intensity to 1 at its maximum. Use complex arithmetic and numpy arrays, avoid Python loops during integration acress the slit. Plot the result. Details of the calculation: Wavelength of the red laser is lambd = 0.6e-6 m. [Don't name the variable lambda, this is a reserved word in Python, and the interpreter won't like it.] Create arrays of 1001 values of x (coordinate across the slit) and y (screen coordinate). Use one loop to step through all y-values. For every y-value, form arrays with integrand and sum them using np.sum() method: Formula for the intensity of light I(y) on the screen: I(y) = |f(y)|**2 where f(y) is the complex-valued amplitude of light, an integral you can perform with any numerical integration scheme you like f(y) = (z**2+y**2)**(-1/2) * integral_{-X}^{+X} exp[i*k*delta] dx, delta = sqrt(z**2+(y-x)**2) - z k = 2 pi / lambd k is the wavenumber, X = 15e-6 m is the half-width of the slit, z = 1. m, y is from -0.06 to +0.06 m. Comment: Delta is the difference of optical paths and its multiplication by wavenumber k=2pi/lambd establishes phase of the light wave on the screen. The complex arithmetic takes care of addition of contributing waves with the correct phase and amplitude. Problem 2 [8pt] - Least Sqares fit to aircraft test data _________________________________________________________ Lectures and textbook dealt with the Least Squares method of fitting theory to data. Here is your chance to use it on aeronautical data. A prototype of a medium-range jetliner was tested in flight, and the data gathered in file http://planets.utsc.utoronto.ca/~pawel/pyth/drag_test_data.dat were obtained. The file looks like this: 7.500e+01 7.136e+04 7.900e+01 7.011e+04 (...) 1.510e+02 9.940e+04 1.600e+02 1.074e+05 1.630e+02 1.109e+05 It contains two columns, one is the measured speed V [m/s], the so-called indicated or calibrated airspeed (in m/s), and the second is the force F [N] of total aerodynamic drag in non-accelerated flight. If you wonder how one can measure such a force, there are a few methods, but the simplest is to let the engines idle, and the plane steadily glide. Force of drag, measured in newtons, can be estimated from the sink rate of the plane. We will return to the minimum sink rate at the end. Read in the data, and fit a two-component theoretical drag curve consisting of a sum of lift-induced drag (https://en.wikipedia.org/wiki/Lift-induced_drag) and parasitic drag (also described/linked on that page). Understanding of all the information on the wiki page is not necessary to solve this problem, but I just thought you might be curious why induced drag falls with the square of the speed v, while the normal (parasitic) aerodynamical drag grows with the square of v. We can write the theoretical dependence as: F(v) = C f0(v) + Ap f1(v) = C/v**2 + Ap (rho v**2/2), (1) where F(v) = total drag force, f_0(v)=v**(-2), and f_1(v)=(rho/2)*v**2, are two known functions of v, rho = sea-level density of air in standard conditions = 1.225 kg/m^3 (this and all other units are S.I.) C = parameter you need to find by fitting data. It is proportional to wing loading factor squared, i.e. (weight of the plane/wingspan)**2, but that is not importan here, let C stand for all unknown factors multiplying v^(-2). C is in S.I. units such that C/v^2 is in newtons. Ap = parameter you need to find by fitting data. This has a meaning of the area of equivalent square plate perpendicular to airflow, which produces the same parasitic drag, so p is a subscript reminding you of that (and of perpendicular plate). Ap is in units of m^2. * * * Least squares method requires solving a matrix equation (here, 2x2 matrix equation) A p = b, where A_ij = sum over all data points (k=0,1,2,..) of f_i(v_k)*f_j(v_k), where i,j={0,1}, v_k is the k-th exprimental value of speed. Vector b is defined as b_i = sum over all data points (k=0,1,2,..) of f_i(v_k)*F_data(v_k). Assemble the matrix A in your code. Also construct the right-hand-side vector b of length 2. Solve the A p = b equation for the vector of parameters p=[C,Ap], which provide the best fit (minimize sum of squares of deviations of theory F(v_k) from observations Fdata(v_k)). You can either solve the simple 2x2 system by hand, or by the np.linalg.solve() method, similarly to how galaxy rotation curves were analyzed in the lecture. The only thing you are not allowed to do is to mindlessly plug the data straight into some function from somebody else's Python module or program, which does the fitting behind the scenes. Once you get the best-fit parameters, use them and create a figure, in which both the experimantal data, the individual components are shown (induced and parasitic drag functions, forming summands in eq. 1), together with the total F(v) fit to data. Prepare the theoretical curve between v=66 m/s and 170 m/s. Although calculations are done in SI units, the values of all forces in N are so large that I recommend feeding them to plt.plot function after division by 9.81e3 (9.81e3 N/T to be precise), thereby converting them from newtons to Tons of force. * * * Finally, analytically (without performing a numerical search), using the found C and Ap values answer the following questions: (i) What is the formula for and the value of v=v_power, the speed at which the plane uses the minimu power in cruise flight. Power = v*F(v). This speed results in minimum fuel flow rate and fuel cost per hour of flight. Also it happens to produce the slowes vertical descent (sink rate) while gliding without engine power. Why? Because sinking is then providing the power. (ii) What is the formula for and the value of v_drag, the speed of minimum total drag force. In other words, F(v_drag) = minimum. At v=v_drag the plane uses the least amount of fuel on a given trajectory leg between points A and B. This is because total energy (total fuel) spent on fighting the dair drag is the drag force times a constant distance from A to B. (iii) Derive the formula for so-called Carson speed of the plane in question, the optimal speed with respect to both time of flight and the fuel used. Since time T is distance/v, and fuel amount scales like F(v), v_T_fuel minimizes the ratio F(v)/v. You need to use calculus on the known F(v) multipied by some power of v, to get all those speeds. Express the 3 optimal speeds in both m/s and km/h, which are more intuitive. * * * Appendix on speed definitions in aeronautics: For your information, here are some clarifications that do not affect your solution, so reading this is optional. The speed we discussed is the IAS or CAS (indicated or calibrated airspeed), not the TAS (true airspeed), which would be closer to actual groundpeed of the airliner. Airplanes travel at more than 160 m/s = 3.6*160 km/h = 576 km/h, ending value in our data. We simply use a different definition of speed, in fact more useful in piloting practice. Here is the difference. The ram pressure (dynamic pressure) p_dyn = (rho_actual/2)*TAS**2 is identical to p_dyn = (rho_std/2)*CAS**2, and both are measured in the plane by the same dynamic pressure gauge, called Pitot tube. Airplane devices are calibrated as if the air was always in standard conditions, and pilots see IAS/CAS not TAS, unless electronically recalculated (air traffic control also uses CAS not TAS). In cruise, rho_actual could be a bit more than 1/3 of the standard density (rho_std) at the see level. Therefore, in the rarified air at cruise altitude of 10-11 km, the plane may move much faster, while working against the same drag force and supporting the same weight: TAS/CAS = (rho_std/rho_actual)**0.5 ~ 3**0.5 ~ 1.7 (For more precise values, use http://www.pdas.com/atmosTable2SI.html). So that's why airliners fly so high: they can move faster, polute less, and provide a less bumpy ride, flying as they are over the clouds. Problem 3 [6pt] - Runge-Kutta integration of Lorenz attractor _____________________________________________________________ Use Runge-Kutta 4th order integration method (programmed by you, do not use any black-box routines) to simulate the dynamical system described by the famous Lorenz set of three nonlinear differential equations (https://en.wikipedia.org/wiki/Lorenz_system) dX/dt = s*(Y-X) dY/dt = X*(r-Z) -Y dZ/dt = X*Y - b*Z Use b = 8/3, s = 10, and r = 28. Start the evolution at t=0 from X,Y,Z = 0, 1, 1.05, and use timestep 0.01 to cover time interval from t=0 to final time T=150. Save all X,Y, and Z values in numpy array as you integrate, for later plotting. After integration completes, plot EITHER 2-d projections onto XY, and XZ plane, OR try to make this prescription for 3-D visualization work for you (either way, you will get full credit if done correctly). Import this module... don't worry, it will write some unimportant warnings: from mpl_toolkits.mplot3d import Axes3D # noqa: F401 unused import # Plot fig = plt.figure(figsize=(12,10)) ax = fig.gca(projection='3d') # now use the values stored in 1-D arrays xs, ys, zs ax.plot(xs, ys, zs, lw=0.5,alpha=0.7) ax.set_xlabel("X Axis") ax.set_ylabel("Y Axis") ax.set_zlabel("Z Axis") ax.set_title("Lorenz Attractor") Whichever way you do the visualization, try to plot an easily visible marker at the end of the curve(s). Then change the initial position a little bit, and observe if the intefgration finished near the previous ending point. How small a variation of initial x, y or z suffices to disturb the ending point unpredictably? This is the famous Butterfly Effect... NOTE: if for some reason you cannot make the RK4 method work, you can perform integration with Euler scheme, for -30% of points. Programs for 2nd order schemes can be found on the web, also we post one on our coding page, but do not use them (do not plagiarize). It's a serious offence. Please read the note at the top of this page about the software that we occasionally use to spot plagiarism. Instead, you must use what hundreds of good companies are using in order to avoid possible litigation from competitors and predecessors. They ask somebody to create a full specification of how the software schould function (get the statement of the problem) and forbidding their engineers to loo at anything else send them to create an independent implementation. They call it "clean room philosphy", and it helped some companies to void major legal troubles. Whole operating systems have been reimplemented cleanly. Problem 4 [5pt] - Uncertainty of parameters of a fit ____________________________________________________ You have some empirical data in file http://planets.utsc.utoronto.ca/~pawel/pyth/hist_data2.dat containig two columns that may be read using np.loadtxt function with delimiter="," option). Ignore the first column which originally used calendar dates. Assume that the second column gives 157 observed y-values y_obs(x), where x are values uniformly distributed between x=0 and x=10. Create the x array. As you know, and have practiced, you can fit this empirical data with a parametrized function, for instance a linear function y(x) = a + b x. If this is a physical science investigation, an important question is not only what values a and b fit the data best, but also what is the uncertainty interval on each parameter, corresponding to the 1-sigma level of uncertainty in the input data. For this purpose, one can write down the solution and use calculus to figure out the sensitivity of a and b to the pertubation of data by +-sigma, a number that is either known from the observation/experiment or can be estimated from deviations between theory and experiment. Alternatively, and somewhat more generally, one can simply create M sets of input data, where each point is subject to a random change with standard deviation = sigma. Notice that +-sigma should be the spread of data points around the inclined line, not around a horizontal line giving the arithmetic mean of all data! (That latter number could be very large, for instance when linear trend is steeply rising, but it has nothing to do with what we care about, which is the mean amount of deviation of data from its general trend, i.e the random error shownig in the data.) After looking at M=10...50 such sets of derived (a,b), you should be able to estimate the uncertainty of a and b reliably enough. Try it! Find the best a and b parameters fitting the data set, plot the best fit, and then estimate the 1-sigma uncertainties in these values by Mte Carlo perturbation of the input data done M=25 times. Explain how much you are purturbing the data and how that corresponds to uncertainly (random error) in original measurements.