Math 171 - Alexiades
                                    Lab 8
                            Finding FWHM

A. Reading data from a file for plotting and exploration :
  Given a (plain text!) file containing N columns of values, we need to read the values into Matlab for exploration and plotting.
  To do this, we need to know the name of the file, and how data is arranged in it, so need to view it in (a plain text!) editor.
1. Write a script "readFile.m" which:
  a. reads a given file "file.dat" into a matrix "TBL",
  b. extracts the columns of the data from "TBL" into arrays "c1", "c2", ..., "cN",
  c. plots some column versus another column, with a title, and labeled axes.
  To test and debug your code, create a "file.dat" with 3 rows and N=2 columns ( so a 3×2 matrix ).
  The first column should be sorted. You can just enter some numbers in "file.dat" in an editor.
  Next, test it on a 4×3 array to make sure it also works. Insert brief comments in you code to document what is being done.

2. To apply this on real data, I am providing a data file, produced by a code that implements a mathematical model for
bacterial chemotaxis (by numerical solution of a system of Partial Differential Equations).
  The file prof300s.dat contains profiles (values at a grid of locations x(i) at some particular time, here at 300sec)
  of several quantities, labeled by column at the top ( x(i)   C(i)   B(i)  etc ). The last column is the index "i" of x(i).
  Of main interest is the 3rd column, labeled "B(i)". It represents bacterial density at location x(i), at the output time
  (normalized, B(i) means B(i)×109 cells/mL). It has a "bump", and want to find its FWHM,
  i.e. the width of the interval [xLeft, xRight] where a horizontal line at half-maximum height intersects the bump.
  Biologically, this is the width of the "band" the bacteria form during chemotaxis; it is visible under a microscope.

B. Plotting  
1. Download the file prof300s.dat and save it as "prof300s.dat" . View it in an editor to see its contents.
  The first two lines contain labels, need to be read over with fgetl .
  The rest of the lines are numerical values to be scanned into a matrix, say "TBL".
  You can see what the B-profile looks like easily in gnuplot:   gnuplot>   plot "prof300s.dat" u 1:3 w lines lw 3

2. Copy your "readFile.m" to "Lab8.m" and modify appropriately to achieve the following:
  a. Read the data from "prof300s.dat" into a matrix "TBL".
  b. Extract each column (which is now a row of "TBL") into a vector (with sensible names: x , C , B , ... ).
    Only need x and B, so may extract only those. Suppress output with ";", they are long!
  c. Find Bmax = maximum value of B, and at which index it occurs, iBmax. Set cut = Bmax/2 .
    Print Bmax and cut nicely (with fprintf). Run your code to make sure it produces values of Bmax, cut.
  d. Plot B vs x, as a red curve. Label the axes. The units of x are micrometers (μm). B is dimensionless.
    Run your code to make sure it produces the plot.
  e. To see where we are cutting the bump, also plot the line at height cut, on the same plot, as a black horizontal line.
    [Note: Once you know the value of "cut", say 1.5, it's easy to do it in Gnuplot (outside of Matlab):
          gnuplot>  cut=1.5;  plot "prof300s.dat" u 1:3 w lines lw 2 , cut lw 3   ]
    To do this in Matlab, will need to create an array, call it "cutArray", of same length as x, with constant values = cut.
    Again, run your code to make sure it produces a plot showing both the bump and the cut-line.
  f. Is this value of cut a good choice here? does it really represent half-maximum of the bump ?
    (note that the "base" value of B is 1, not zero). Figure out a better expression for setting cut.

C. Find FWHM numerically, from arrays of points.
1. Create a Matlab function in a file "fwhm.m", which accepts two arrays (x, y) as input, and returns FWHM of y,
  by implementing the following algorithm (discussed in Part 5 of HW1).
  • Find the maximum value Peak of y, and the index where it occurs, iPeak .
  • Set the cutoff level at half maximum:  cut = Peak/2 , or more appropriately from B.2.f.
  • Find the first point xi (before iPeak) at which yi ≥ cut. Call it xLeft.
  • Find the first xi (after iPeak) at which yi ≤ cut. Call it xRight.
  • Then FWHM = xRight − xLeft .

    2. Modify your "Lab8.m" code to:
      a. Call the fwhm function with arguments (x, B) to find the FWHM of the B array.
      b. Then print it out (nicely, with fprintf) to say:  
        "FWHM of bacteria at time 300 sec is"  FWHM_value  "microns, at height"  cut_value

    D. Prepare Lab8.m for submission:
  • As always, at the top of your Lab8.m code, insert % Name, Date, Lab8
  • At the end of the code insert a separator:   %=====================================%
  • Insert the "fwhm.m" code into "Lab8.m" as a function subprogram, and then another separator.
  • Clean it up. Suppress output with ";".   Insert appropriate comments.
  • Make sure your code runs. It should produce only the plot and FWHM.

  • Publish it in Matlab (works in Matlab, might also work in Octave): >> publish( 'Lab8' , 'pdf' )
      This runs the code and records both the code and its output (thus also the plot) into a file Lab8.pdf
      It will be in the html subfolder. View the PDF to make sure it contains code and plot and FWHM.

    Then submit Lab8.pdf on Canvas by due date.