(************** Content-type: application/mathematica **************
                     CreatedBy='Mathematica 5.0'

                    Mathematica-Compatible Notebook

This notebook can be used with any Mathematica-compatible
application, such as Mathematica, MathReader or Publicon. The data
for the notebook starts with the line containing stars above.

To get the notebook into a Mathematica-compatible application, do
one of the following:

* Save the data starting with the line of stars above into a file
  with a name ending in .nb, then open the file inside the
  application;

* Copy the data starting with the line of stars above to the
  clipboard, then use the Paste menu command inside the application.

Data for notebooks contains only printable 7-bit ASCII and can be
sent directly in email or through ftp in text mode.  Newlines can be
CR, LF or CRLF (Unix, Macintosh or MS-DOS style).

NOTE: If you modify the data for this notebook not in a Mathematica-
compatible application, you must delete the line below containing
the word CacheID, otherwise Mathematica-compatible applications may
try to use invalid cache data.

For more information on notebooks and Mathematica-compatible 
applications, contact Wolfram Research:
  web: http://www.wolfram.com
  email: info@wolfram.com
  phone: +1-217-398-0700 (U.S.)

Notebook reader applications are available free of charge from 
Wolfram Research.
*******************************************************************)

(*CacheID: 232*)


(*NotebookFileLineBreakTest
NotebookFileLineBreakTest*)
(*NotebookOptionsPosition[     89499,       1802]*)
(*NotebookOutlinePosition[     90145,       1824]*)
(*  CellTagsIndexPosition[     90101,       1820]*)
(*WindowFrame->Normal*)



Notebook[{
Cell[BoxData[
    \(\(\(\ \)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\( (*\ 
      Routine\ for\ Data\ Reduction\ of\ Seeing\ data\  - \ 
        on\ the\ wave\ crest\ *) \)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \( (*\ 
      The\ first\ part\ of\ the\ data\ analysis\ will\ look\ at\ the\ flat - 
        calibrated\ raw\ data\ coming\ from\ the\ CCD\ \((input\ of\ text\ \
file\ in\ ReadList)\); \ 
      the\ idea\ is\ to\ have\ a\ good\ 3  D\ view\ of\ it, \ 
      where\ the\ vertical\ axis\ are\ the\ ADUs\ of\ each\ pixel\ *) \)], \
"Input"],

Cell[BoxData[
    \(\(\( (*\ 1  st, \ 
      we\ need\ to\ input\ the\ matrix\ \((x, 
          y)\)\ size\ of\ the\ CCD\ chip\ *) \)\(\[IndentingNewLine]\)\( (*\ 
      Here, \ and\ elsewhere\ in\ the\ routine, \ 
      relevant\ numbers\ will\ be\ printed\ our\ so\ that\ after\ saving\ the\
\ file\ these\ can\ be\ looked\ up\ later; \ 
      we\ start\ doing\ so\ for\ the\ [x, 
          y]\ CCD\ chip\ size\ *) \)\)\)], "Input"],

Cell[BoxData[
    \(matx = 
      Input["\<Please enter the x-direction number of pixels in the CCD chip:\
\>"]; maty = 
      Input["\<Please enter the y-direction number of pixels in the CCD chip:\
\>"]; Print["\<CCD chip size (x,y): (\>", matx, "\<,\>", 
      maty, "\<)\>"];\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ Routine\ 1\  - \ 
        Initial\ look\ at\ the\ data\ *) \)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \( (*\ The\ 1  
        st\ plot\ will\ give\ an\ idea\ of\ where\ the\ data\ lie\ in\ the\ \
\((x, y, ADU)\)\ space; \ 
      note\ that\ the\ initial\ ADU\ maximum\ is\ set\ at\ 60000, \ 
      but\ this\ might\ be\ changed\ higher\ \((if\ the\ data\ come\ out\ \
chopped)\)\ or\ lower\ \((closer\ to\ the\ maximum\ ADU\ value)\); \ 
      the\ maximum\ value\ of\ ADUs, \ 
      that\ the\ routine\ also\ gives\ as\ output, \ 
      helps\ in\ this\ matter; \ 
      only\ change\ the\ viewpoint\ as\ a\ last\ resource\ *) \)], "Input"],

Cell[BoxData[
    \( (*\ The\ 1  
          st\ instruction\ is\ actually\ to\ ignore\ spell\ checking\ since\ \
it\ can\ be\ annoying\ for\ similarly - named\ variables\ *) \)], "Input"],

Cell[BoxData[
    StyleBox[\(Off[General::spell1]\),
      "MR"]], "Input"],

Cell[BoxData[
    \( (*\ Obviously, \ 
      change\ the\ SetDirectory\ according\ to\ the\ place\ where\ you\ have\ \
the\ CCD\ text\ file; \ 
      and\ place\ its\ name\ on\ "\<ReadList\>"\ *) \)], "Input"],

Cell[BoxData[{
    \(SetDirectory["\<C:\Users\FullDirectoryPathHERE\>"]; \ 
    a = ReadList["\<FILE_NAME_HERE.txt\>", 
        Table[Number, {matx}]];\), "\[IndentingNewLine]", 
    \(ListPlot3D[a, ColorFunction \[Rule] Hue, 
      PlotRange -> {{0, matx}, {0, maty}, {0, 60000}}, 
      ViewPoint \[Rule] {1, 1,  .2}]; m = Max[a]; 
    Print["\<Peak value (ADU):  \>", m]\)}], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ Routine\ 2\  - \ 
        Interact\ with\ the\ \(data : \ 
            final\ graph\ production\ and\ view\)\ *) \)\(\[IndentingNewLine]\
\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \( (*\ Given\ the\ peak\ value, \ 
      the\ following\ routine\ will\ have\ to\ be\ adjusted\ until\ we\ have\ \
a\ "\<nice looking\>"\ 3  D\ plot\ of\ the\ Polaris\ path; \ 
      we\ usually\ limit\ the\ maximum\ to\ the\ next\ "\<ten of thousand\>"\ \
and\ guess\ a\ tighter\ SQUARE\ \((x1, x2)\) - \((y1, 
            y2)\)\ matrix\ from\ the\ plot\ above, \ 
      eventually\ enlarging\ \(it : \ trial\ an\ error\); \ in\ the\ end, \ 
      changing\ the\ viewpoint\ might\ be\ necessary\ too\ *) \)], "Input"],

Cell[BoxData[
    \( (*\ Once\ you\ have\ settled\ on\ a\ fine\ looking\ plot, \ 
      the\ final\ values\ must\ be\ inserted\ by\ running\ the\ routine\ that\
\ follows; \ the\ print\ out\ confirms\ them\ *) \)], "Input"],

Cell[BoxData[
    \(nx1 = Input["\<x1?\>"]; nx2 = Input["\<x2?\>"]; ny1 = Input["\<y1?\>"]; 
    ny2 = Input["\<y2?\>"]; 
    Print["\<x1= \>", nx1, "\<, x2= \>", nx2, "\<, y1= \>", \ 
      ny1, \ "\<, y2= \>", \ ny2];\)], "Input"],

Cell[BoxData[
    \( (*\ The\ final\ \((x1 - x2)\) - \((y1 - 
              y2)\)\ values\ are\ now\ inserted, \ 
      and\ there\ we\ \(go!\)\ \ *) \)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ Routine\ 3\ \((Main)\)\  - \ 
        Organizing\ the\ data\ into\ 2  x2\ pixels, \ 
      preparing\ vectors\ c\ and\ f\ for\ analysis; \ be\ patient, \ 
      it\ might\ take\ a\ few\ \(minutes!\)\ *) \)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \( (*\ The\ following\ routine\ will\ identify\ the\ coordinates\ \((x, 
          y)\)\ of\ the\ maximum\ value\ in\ ADU\ \((m)\); \ 
      right\ now\ we\ also\ define, \ from\ this\ value, \ 
      the\ maximum\ range\ for\ the\ contour\ plots\  - \ 
        topgr\ *) \)], "Input"],

Cell[BoxData[
    \(b = Transpose[a]; \ \ x = 1; y = 1; mx = 0; my = 0; 
    While[y < maty + 1, \ 
      While[x < matx + 1, \ 
        If[b[\([x, y]\)] == m, mx = x; my = y; x = matx + 1; y = maty + 1, 
          x = x + 1, Print["\<erro\>"]]]; x = 1; y = y + 1]; \ 
    topgr = m + 5000;\)], "Input"],

Cell[BoxData[
    \( (*\ Now, \ starting\ with\ the\ maximum\ value, \ 
      we\ define\ the\ 1  st\ 2  x2\ ' group'\ which\ has\ the\ 1  
        st\ maximum, \ 2  nd\ maximum\ \((close)\), \ 3  rd\ maximum, \ 
      etc . \ all\ together\  - \ 
        other\ values\ further\ away\ may\ be\ higher\ but\ we\ use\ the\ \
peak\ as\ an\ ' anchor'\ for\ defining\ the\ "\<peak 4x4 group\>"; \ 
      the\ 2  x2\ peak\ will\ be\ defined\ within\ the\ eight\ pixels\ that\ \
surround\ the\ peak\ *) \)], "Input"],

Cell[BoxData[
    \( (*\ We\ first\ find\ the\ 2  nd\ maximum - m2\ *) \)], "Input"],

Cell[BoxData[
    \(x = mx - 1; y = my - 1; i = 0; j = 0; u = 100; 
    While[y < my + 2, 
      While[x < mx + 2, 
        If[b[\([x, y]\)] > u, u = b[\([x, y]\)]; i = x; j = y]; x = x + 1]; 
      y = y + 2; x = mx - 1]; 
    If[b[\([mx - 1, my]\)] > u, u = b[\([mx - 1, my]\)]; i = mx - 1; 
      j = my]; \ 
    If[b[\([mx + 1, my]\)] > u, u = b[\([mx + 1, my]\)]; i = mx + 1; j = my]; 
    m2 = u; \ m2x = i; m2y = j;\)], "Input"],

Cell[BoxData[
    \(\(\( (*\ Now, \ 
      if\ the\ 2  nd\ maximum\ is\ diagonal\ as\ regards\ the\ peak\ m, \ 
      the\ 2  x2\ peak\ is\ immediately\ defined; \ 
      we\ check\ this\ next\ ; \ 
      note\ that\ we\ define\ the\ 2  x2\ peak\ as\ a\ new\ pixel\ in\ the\ 1 
         st\ element\ of\ a\ new\ matrix\ d\ where\ the\ ADU\ is\ the\ \
average\ of\ the\ four\ values\ and\ the\ coordinates\ of\ the\ 2  
        x2\ pixel\ are\ the\ ones\ of\ the\ point\ at\ the\ intersection\ of\ \
the\ four\ pixels; \ similarly, \ 
      another\ matrix\ e\ is\ built\ \((1  
            st\ element)\)\ containing\ the\ standard\ deviations\ of\ the\ \
calculus\ *) \)\(\ \)\( (*\ 
      Now\ inserted\ into\ the\ 3  
        rd\ IF\ of\ the\ one\ above\ \((option\ if\ test\ is\ false)\)\ we\ \
work\ the\ more\ difficult\ situation\ of\ the\ 2  
        nd\ maximum\ not\ being\ in\ a\ diagonal; \ in\ this\ case, \ 
      it\ is\ either\ left, \ up, \ right\ or\ down\ from\ the\ main\ peak; \ 
      to\ decide, \ 
      we\ will\ have\ to\ look\ for\ the\ 3  
        rd\ peak\ in\ the\ other\ two\ pixels\ that\ will\ complete\ the\ 2  
        x2; \ two\ hypothesis\ in\ each\ case*) \)\)\)], "Input"],

Cell[BoxData[
    \(esq = 0; drt = 0; cim = 0; bxo = 0; dx = 0; dy = 0; dm = 0; es = 0; 
    d = {}; e = {}; i = 0; j = 0; 
    If[Abs[mx - m2x] \[Equal] 1, i = 1, i = 0]; \ 
    If[Abs[my - m2y] \[Equal] 1, j = 1, j = 0]; 
    If[i + j \[Equal] 2, 
      If[mx > m2x && m2y > my, dx = mx - 0.5; dy = my + 0.5; 
        dm = N[Mean[{m, m2, b[\([mx - 1, my]\)], b[\([mx, my + 1]\)]}]]; 
        es = N[StandardDeviation[{m, m2, b[\([mx - 1, my]\)], 
                b[\([mx, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
        e = Append[e, {dx, dy, es}], 
        If[m2x > mx && m2y > my, dx = mx + 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my]\)], b[\([mx, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], 
          If[m2x > mx && my > m2y, dx = mx + 0.5; dy = my - 0.5; 
            dm = N[Mean[{m, m2, b[\([mx + 1, my]\)], b[\([mx, my - 1]\)]}]]; 
            es = N[StandardDeviation[{m, m2, b[\([mx + 1, my]\)], 
                    b[\([mx, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
            e = Append[e, {dx, dy, es}], 
            If[mx > m2x && my > m2y, dx = mx - 0.5; dy = my - 0.5; 
              dm = N[Mean[{m, m2, b[\([mx - 1, my]\)], 
                      b[\([mx, my - 1]\)]}]]; 
              es = N[StandardDeviation[{m, m2, b[\([mx - 1, my]\)], 
                      b[\([mx, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
              e = Append[e, {dx, dy, es}]]]]], 
      If[m2x == mx && m2y > my, 
        esq = Max[b[\([mx - 1, my + 1]\)], b[\([mx - 1, my]\)]]; 
        dir = Max[b[\([mx + 1, my + 1]\)], b[\([mx + 1, my]\)]]; 
        If[esq > dir, dx = mx - 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx - 1, my]\)], 
                  b[\([mx - 1, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx - 1, my]\)], 
                  b[\([mx - 1, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], dx = mx + 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx + 1, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx + 1, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}]]]; 
      If[m2x == mx && m2y < my, 
        esq = Max[b[\([mx - 1, my - 1]\)], b[\([mx - 1, my]\)]]; 
        dir = Max[b[\([mx + 1, my - 1]\)], b[\([mx + 1, my]\)]]; 
        If[esq > dir, dx = mx - 0.5; dy = my - 0.5; 
          dm = N[Mean[{m, m2, b[\([mx - 1, my]\)], 
                  b[\([mx - 1, my - 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx - 1, my]\)], 
                  b[\([mx - 1, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], dx = mx + 0.5; dy = my - 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx + 1, my - 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx + 1, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}]]]; 
      If[m2y == my && m2x < mx, 
        cim = Max[b[\([mx - 1, my + 1]\)], b[\([mx, my + 1]\)]]; 
        bxo = Max[b[\([mx - 1, my - 1]\)], b[\([mx, my - 1]\)]]; 
        If[cim > bxo, dx = mx - 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx - 1, my + 1]\)], 
                  b[\([mx, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx - 1, my + 1]\)], 
                  b[\([mx, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], dx = mx - 0.5; dy = my - 0.5; 
          dm = N[Mean[{m, m2, b[\([mx - 1, my - 1]\)], 
                  b[\([mx, my - 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx - 1, my - 1]\)], 
                  b[\([mx, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}]]]; 
      If[m2y == my && m2x > mx, 
        cim = Max[b[\([mx + 1, my + 1]\)], b[\([mx, my + 1]\)]]; 
        bxo = Max[b[\([mx + 1, my - 1]\)], b[\([mx, my - 1]\)]]; 
        If[cim > bxo, dx = mx + 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my + 1]\)], 
                  b[\([mx, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my + 1]\)], 
                  b[\([mx, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], dx = mx + 0.5; dy = my - 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my - 1]\)], 
                  b[\([mx, my - 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my - 1]\)], 
                  b[\([mx, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}]]]];\)], "Input"],

Cell[BoxData[
    \(\(\( (*\ Now\ we\ organize\ all\ data\ into\ 2  x2\ pixels; \ 
      special\ care\ has\ to\ be\ taken\ when\ reaching\ the\ corners\ and\ \
this\ is\ why\ the\ previous\ work\ was\ \(done : \ 
          peak\ identification\ allows\ to\ decide\ where\ to\ chop\ the\ \
data\ in\ x\ and\ y\ in\ order\ to\ ensure\ that\ both\ axis\ have\ a\ number\
\ of\ pixels\ divisible\ by\ four\); \ 
      note\ that\ we\ initialize\ the\ matrices\ d\ and\ e\ so\ that\ they\ \
are\ ready\ to\ remap\ the\ CCD\ chip\ into\ 2  x2\ ; \ 
      we\ also\ create\ vector\ d1\ \((and\ its\ transpose\ d3)\)\ so\ that\ \
we\ can\ later\ plot\ the\ fit\ on\ top\ of\ the\ 2  x2\ data, \ 
      as\ well\ as\ on\ top\ of\ the\ original\ 1  x1\ data*) \)\(\ \)\( (*\ 
      this\ routine\ takes\ a\ few\ minutes\ to\ run\ *) \)\)\)], "Input"],

Cell[BoxData[
    \(mx = dx; my = dy; m = dm; ordd = 0; ul = 0; ur = 0; dl = 0; dr = 0; 
    i = dx - 0.5; j = dy - 0.5; dx = 0; dy = 0; dm = 0; es = 0; xlo = 1; 
    xhi = matx; ylo = 1; yhi = maty; d = {}; e = {}; d1 = {}; \ d2 = {}; \ 
    If[EvenQ[IntegerPart[j]] \[Equal] True, ylo = 2; yhi = maty - 1]; 
    If[EvenQ[IntegerPart[i]] \[Equal] False, xhi = matx - 1, xlo = 2]; 
    y = ylo; v = yhi; x = xlo; l = xhi; 
    While[y < v + 1, 
      While[x < l + 1, dl = b[\([x, y]\)]; dr = b[\([x + 1, y]\)]; 
        ul = b[\([x, y + 1]\)]; ur = b[\([x + 1, y + 1]\)]; dx = x + 0.5; 
        dy = y + 0.5; dm = N[Mean[{dl, dr, ul, ur}]]; 
        es = N[StandardDeviation[{dl, dr, ul, ur}]]; 
        d = Append[d, {dx, dy, dm}]; d2 = Append[d2, dm]; 
        d2 = Append[d2, dm]; 
        If[dx == mx && dy \[Equal] my, ordd = Length[d]]; 
        e = Append[e, {dx, dy, es}]; x = x + 2]; x = xlo; 
      d1 = Append[d1, d2]; d1 = Append[d1, d2]; d2 = {}; y = y + 2]; 
    d3 = Transpose[d1];\)], "Input"],

Cell[BoxData[
    \( (*\ Now, \ 
      we\ build\ the\ main\ \(routine : \ 
          starting\ from\ the\ peak\ and\ building\ up\ the\ path\ both\ ways\
\), \ until\ reaching\ a\ value\  < 10000\ ADU; \ 
      always\ picking\ the\ highest - valued\ pixel\ around*) \)], "Input"],

Cell[BoxData[
    \( (*\ We\ first\ do\ it\ for\ the\ peak\ only; \ 
      this\ will\ define\ the\ direction\ of\ the\ initial\ path\ \((2  
              nd\ maximum - m2)\); \ 
      the\ second\ similar\ subroutine\ does\ the\ other\ bit\ \((3  
              rd\ maximum - 
            m3)\)\ that\ will\ be\ used\ to\ complete\ the\ path\ later, \ 
      in\ the\ opposite\ direction\ *) \)], "Input"],

Cell[BoxData[
    \(c = {{mx, my, m}}; k = 0; l = 0; m2xyADU = 0; u = 10001; 
    ordem = ordd - IntegerPart[\((matx + 3)\)/2]; 
    While[l < 2, 
      While[k < 3, ordem = ordem + 1; 
        If[d[\([ordem, 3]\)] > 10000 && d[\([ordem, 3]\)] > u, 
          u = d[\([ordem, 3]\)]; m2xyADU = ordem]; k = k + 1]; k = 0; 
      ordem = ordem + matx - 4; l = l + 1]; ordem = ordd - 3; 
    While[k < 2, ordem = ordem + 2; 
      If[d[\([ordem, 3]\)] > 10000 && d[\([ordem, 3]\)] > u, 
        u = d[\([ordem, 3]\)]; m2xyADU = ordem]; k = k + 1]; m2 = u; 
    m2x = d[\([m2xyADU, 1]\)]; m2y = d[\([m2xyADU, 2]\)]; 
    If[m2 > 10000, c = Append[c, {m2x, m2y, m2}]];\)], "Input"],

Cell[BoxData[
    \(k = 0; l = 0; m3xyADU = 0; u = 10001; 
    ordem = ordd - IntegerPart[\((matx + 3)\)/2]; f = {}; 
    While[l < 2, 
      While[k < 3, ordem = ordem + 1; 
        If[d[\([ordem, 3]\)] > 10000 && d[\([ordem, 3]\)] > u && 
            d[\([ordem, 1]\)] \[NotEqual] m2x && 
            d[\([ordem, 2]\)] \[NotEqual] m2y, u = d[\([ordem, 3]\)]; 
          m3xyADU = ordem]; k = k + 1]; k = 0; ordem = ordem + matx - 4; 
      l = l + 1]; ordem = ordd - 3; 
    While[k < 2, ordem = ordem + 2; 
      If[d[\([ordem, 3]\)] > 10000 && d[\([ordem, 3]\)] > u && 
          d[\([ordem, 1]\)] \[NotEqual] m2x && 
          d[\([ordem, 2]\)] \[NotEqual] m2y, u = d[\([ordem, 3]\)]; 
        m3xyADU = ordem]; k = k + 1]; m3 = u; m3x = d[\([m3xyADU, 1]\)]; 
    m3y = d[\([m3xyADU, 2]\)]; 
    If[m3 > 10000, f = Append[f, {m3x, m3y, m3}]];\)], "Input"],

Cell[BoxData[
    \( (*\ 
      We\ now\ have\ the\ main\ routine\ \((currently\ limited\ to\ 100\ \
cycles, \ 
          otherwise\ it\ might\ enter\ an\ infinite\ cycle\ mode\ before\ \
reaching\ low\ pixels; \ \ \(if\ it\ reaches\ values\  < 
                10000\ ADU\ before\ the\ 100\ cycles\ it\ will\ conclude\ \
with\ \ error \((s)\) ... \)\ [don' t\ worry])\); \ 
      this\ will\ start\ from\ the\ 2  
        nd\ maximum\ and\ explore\ only\ three\ pixels; \ 
      if\ m - m2\ is\ a\ diagonal, \ these\ are\ the\ "\<corner\>"\ arrow; \ 
      otherwise\ they\ will\ be\ in\ a\ line; \ 
      the\ routine\ stops\ when\ reaching\ a\ value\ lower\ than\ 10000\ \
\((st = 1)\); \ 
      then\ we\ will\ do\ the\ same\ in\ a\ separate\ routine\ \((for\ now)\)\
\ for\ the\ third\ maximum\ \((m3)\)\ so\ that\ we\ explore\ the\ other\ \
direction\ from\ the\ peak\ \((m)\); \ 
      a\ later\ graph\ inspection\ might\ justify\ running\ the\ routine\ in\ \
"\<chunks\>"\ \((if\ there\ are\ some\  < 
              10000\ ADU\ areas\ that\ split\ higher\ values\ ones)\) : \ 
        what\ we\ call\ "\<lobe work\>"\ *) \)], "Input"],

Cell[BoxData[{
    \(st = 0; temporario = 0; ciclos1 = 0; m1xyADU = ordd; m2ADU = m2xyADU; 
    val382 = IntegerPart[\((matx - 1)\)/2];\), "\[IndentingNewLine]", 
    \(\(While[st == 0 && ciclos1 < 100, orddif = m2xyADU - ordd; 
        tempor = m2xyADU; u = 0; chave = Abs[val382 - Abs[orddif]]; 
        ordem = 0; teste = 0; 
        If[d[\([ordd + 2*orddif, 3]\)] > 10000 && 
            d[\([ordd + 2*orddif, 3]\)] > u, u = d[\([ordd + 2*orddif, 3]\)]; 
          ordem = ordd + 2*orddif]; \[IndentingNewLine]If\ [chave == 1, 
          teste = orddif - val382; 
          If[teste \[Equal] 1 || teste \[Equal] \(-matx\) + 2, 
            If[d[\([ordd + 2*orddif - 1, 3]\)] > 10000 && 
                d[\([ordd + 2*orddif - 1, 3]\)] > u, 
              u = d[\([ordd + 2*orddif - 1, 3]\)]; 
              ordem = ordd + 2*orddif - 1]; 
            If[d[\([ordd + orddif + 1, 3]\)] > 10000 && 
                d[\([ordd + orddif + 1, 3]\)] > u, 
              u = d[\([ordd + orddif + 1, 3]\)]; ordem = ordd + orddif + 1], 
            If[d[\([ordd + 2*orddif + 1, 3]\)] > 10000 && 
                d[\([ordd + 2*orddif + 1, 3]\)] > u, 
              u = d[\([ordd + 2*orddif + 1, 3]\)]; 
              ordem = ordd + 2*orddif + 1]; 
            If[d[\([ordd + orddif - 1, 3]\)] > 10000 && 
                d[\([ordd + orddif - 1, 3]\)] > u, 
              u = d[\([ordd + orddif - 1, 3]\)]; ordem = ordd + orddif - 1]], 
          If[chave \[Equal] 0, 
            If[d[\([ordd + 2*orddif + 1, 3]\)] > 10000 && 
                d[\([ordd + 2*orddif + 1, 3]\)] > u, 
              u = d[\([ordd + 2*orddif + 1, 3]\)]; 
              ordem = ordd + 2*orddif + 1]; 
            If[d[\([ordd + 2*orddif - 1, 3]\)] > 10000 && 
                d[\([ordd + 2*orddif - 1, 3]\)] > u, 
              u = d[\([ordd + 2*orddif - 1, 3]\)]; 
              ordem = ordd + 2*orddif - 1], 
            If[d[\([tempor + \((val382 + 1)\)*orddif, 3]\)] > 10000 && 
                d[\([tempor + \((val382 + 1)\)*orddif, 3]\)] > u, 
              u = d[\([tempor + \((val382 + 1)\)*orddif, 3]\)]; 
              ordem = tempor + \((val382 + 1)\)*orddif]; 
            If[d[\([tempor - \((val382 - 1)\)*orddif, 3]\)] > 10000 && 
                d[\([tempor - \((val382 - 1)\)*orddif, 3]\)] > u, 
              u = d[\([tempor - \((val382 - 1)\)*orddif, 3]\)]; 
              ordem = tempor - \((val382 - 1)\)*orddif]]]; m2 = u; 
        m2x = d[\([ordem, 1]\)]; m2y = d[\([ordem, 2]\)]; m2xyADU = ordem; 
        ordd = tempor; ciclos1 = ciclos1 + 1; 
        If[m2 > 10000, c = Append[c, {m2x, m2y, m2}], st = 1]];\)\)}], "Input"],

Cell[BoxData[
    \(\(\( (*\ We\ now\ do\ the\ main\ routine\ for\ the\ 3  rd\ maximum, \ 
      exactly\ the\ same\ way\ ; \ however, \ 
      we\ will\ build\ vector\ f\ so\ that\ we\ can\ later\ combine\ with\ \
vector\ c\ to\ make\ a\ single\ path\ for\ the\ Pole\ Star\ from\ our\ data\ \
\((c + f)\)\ *) \)\(\ \)\(\[IndentingNewLine]\)\( (*\ once\ again, \ 
      error \((s)\)\ is\ ok\ *) \)\)\)], "Input"],

Cell[BoxData[{
    \(st = 0; ciclos2 = 0; orddif = 0; temporario = 0; ordd = m1xyADU; 
    m3ADU = m3xyADU;\), "\[IndentingNewLine]", 
    \(\(While[st == 0 && ciclos2 < 100, orddif = m3xyADU - ordd; 
        temporario = m3xyADU; u = 0; chave = Abs[val382 - Abs[orddif]]; 
        ordem = 0; teste = 0; 
        If[d[\([ordd + 2*orddif, 3]\)] > 10000 && 
            d[\([ordd + 2*orddif, 3]\)] > u, u = d[\([ordd + 2*orddif, 3]\)]; 
          ordem = ordd + 2*orddif]; \[IndentingNewLine]If\ [chave == 1, 
          teste = orddif - val382; 
          If[teste \[Equal] 1 || teste \[Equal] \(-matx\) + 2, 
            If[d[\([ordd + 2*orddif - 1, 3]\)] > 10000 && 
                d[\([ordd + 2*orddif - 1, 3]\)] > u, 
              u = d[\([ordd + 2*orddif - 1, 3]\)]; 
              ordem = ordd + 2*orddif - 1]; 
            If[d[\([ordd + orddif + 1, 3]\)] > 10000 && 
                d[\([ordd + orddif + 1, 3]\)] > u, 
              u = d[\([ordd + orddif + 1, 3]\)]; ordem = ordd + orddif + 1], 
            If[d[\([ordd + 2*orddif + 1, 3]\)] > 10000 && 
                d[\([ordd + 2*orddif + 1, 3]\)] > u, 
              u = d[\([ordd + 2*orddif + 1, 3]\)]; 
              ordem = ordd + 2*orddif + 1]; 
            If[d[\([ordd + orddif - 1, 3]\)] > 10000 && 
                d[\([ordd + orddif - 1, 3]\)] > u, 
              u = d[\([ordd + orddif - 1, 3]\)]; ordem = ordd + orddif - 1]], 
          If[chave \[Equal] 0, 
            If[d[\([ordd + 2*orddif + 1, 3]\)] > 10000 && 
                d[\([ordd + 2*orddif + 1, 3]\)] > u, 
              u = d[\([ordd + 2*orddif + 1, 3]\)]; 
              ordem = ordd + 2*orddif + 1]; 
            If[d[\([ordd + 2*orddif - 1, 3]\)] > 10000 && 
                d[\([ordd + 2*orddif - 1, 3]\)] > u, 
              u = d[\([ordd + 2*orddif - 1, 3]\)]; 
              ordem = ordd + 2*orddif - 1], 
            If[d[\([temporario + \((val382 + 1)\)*orddif, 3]\)] > 10000 && 
                d[\([temporario + \((val382 + 1)\)*orddif, 3]\)] > u, 
              u = d[\([temporario + \((val382 + 1)\)*orddif, 3]\)]; 
              ordem = temporario + \((val382 + 1)\)*orddif]; 
            If[d[\([temporario - \((val382 - 1)\)*orddif, 3]\)] > 10000 && 
                d[\([temporario - \((val382 - 1)\)*orddif, 3]\)] > u, 
              u = d[\([temporario - \((val382 - 1)\)*orddif, 3]\)]; 
              ordem = temporario - \((val382 - 1)\)*orddif]]]; m3 = u; 
        m3x = d[\([ordem, 1]\)]; m3y = d[\([ordem, 2]\)]; m3xyADU = ordem; 
        ordd = temporario; ciclos2 = ciclos2 + 1; 
        If[m3 > 10000, f = Append[f, {m3x, m3y, m3}], st = 1]];\)\)}], "Input"],

Cell[BoxData[
    \( (*\ There\ is\ a\ final\ check\ to\ make, \ 
      before\ establishing\ the\ final\ complete\ \(path : \ 
          visually\ look\ \((ideally, \ print)\)\ at\ both\ \((x, 
              y)\)\ graphs\ of\ c\ and\ f\ to\ make\ sure\ they\ will\ \
complete\ each\ other\ logically\); \ 
      we\ may\ decide\ new\ sections\ to\ remove; \ in\ that\ case, \ 
      just\ go\ to\ routine\ 4\ next\ and\ change\ the\ chopping\ pixel; \ 
      for\ reference, \ 
      the\ peak\ and\ 2  
        nd\ peak\ \ have\ coordinates\ printed\ out\ here\ \((1  
            st\ graph)\), \ together\ with\ vector\ c\ in\ full; \ 
      also\ the\ 3  rd\ peak, \ 
      in\ its\ own\ graph\ \((and\ vector\ f\ in\ full)\)\ *) \)], "Input"],

Cell[BoxData[
    \(l = 1; e = {}; n2 = Length[c]; 
    While[l \[LessEqual] n2, e = Append[e, {c[\([l, 1]\)], c[\([l, 2]\)]}]; 
      l = l + 1]; 
    ListPlot[e, PlotRange \[Rule] All, PlotStyle -> PointSize[0.013]]; \ 
    Print["\<Coordinates of the peak (x,y):  \>", 
      d[\([m1xyADU, 1]\)], "\<   \>", d[\([m1xyADU, 2]\)]]; 
    Print["\<Coordinates of the 2nd maximum (x,y):  \>", 
      d[\([m2ADU, 1]\)], "\<   \>", d[\([m2ADU, 2]\)]]; 
    Print["\< Vector C:  \>", c]\)], "Input"],

Cell[BoxData[
    \(l = 1; e = {}; n2 = Length[f]; 
    While[l \[LessEqual] n2, e = Append[e, {f[\([l, 1]\)], f[\([l, 2]\)]}]; 
      l = l + 1]; 
    ListPlot[e, PlotRange \[Rule] All, PlotStyle -> PointSize[0.013]]; 
    Print["\<Coordinates of the 3rd maximum (x,y):  \>", 
      d[\([m3ADU, 1]\)], "\<   \>", d[\([m3ADU, 2]\)]]; 
    Print["\< Vector F:  \>", f]\)], "Input"],

Cell[BoxData[
    \( (*\ Note\ that, \ before\ proceeding, \ 
      clean\ versions\ of\ c\ and\ f\ must\ be\ created\ if\ the\ cycles\ \
reached\ the\ number\ 100 : \ 
        we\ have\ to\ remove\ all\ "\<repeat\>"\ sections\ ; \ \ the\ maximum\
\ point\ of\ interest\ \((before\ repeats)\)\ MUST\ BE\ CHANGED\ IN\ THE\ \
ROUTINE\ FOR\ EACH\ \(CASE!\)\ *) \)], "Input"],

Cell[BoxData[
    \( (*\ So, \ we\ first\ run\ the\ following\ routine, \ 
      to\ check\ if\ that\ was\ the\ case; \ 
      if\ one\ or\ two\ of\ the\ cycles\ reached\ 100\ then\ the\ \
corresponding\ routine\ 4\ that\ follows\ should\ be\ run\ to\ inspect\ the\ \
graph \((s)\)\ of\ the\ relevant\ vector \((s)\)*) \)], "Input"],

Cell[BoxData[
    \(Print["\<Cycles of vector c, which includes the peak:  \>", 
      ciclos1 + 1]; 
    Print["\<Cycles of vector f, which does not include the peak:  \>", 
      ciclos2];\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\( (*\ 
      Note\ that\ in\ the\ example\ used\ here, \ 
      neither\ vector\ reached\ 100\ cycles\ so\ we\ would\ jump\ straight\ \
to\ routine\ 5; \ however, \ 
      routine\ 4\ is\ also\ used\ to\ redefine\ one\ or\ two\ of\ vectors\ c/
          f; \ since, \ in\ this\ example, \ we\ want\ to\ do\ that, \ 
      we\ will\ actually\ use\ it\ \ *) \)\(\[IndentingNewLine]\)\( (*\ 
      In\ this\ example\ we\ will\ remove\ the\ last\ point\ of\ "\<f\>" + 
        c\ , {427.5, 376.5, 10671.5}, \ 
      due\ to\ a\ 59  %\ drop\ *) \)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\( (*\ 
      Routine\ 4\ \((Cycles)\)\  - \ 
        eventually\ redefining\ vector\ c\ and/
            or\ vector\ f\ \
*) \)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\( (*\ 
      This\ routine\ might\ be\ used\ after\ graph\ inspection\ and/
          or\ if\ any\ of\ the\ cycles\ for\ vector\ c/
          f\ reached\ 100\ \
*) \)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \(lma = 
      Input["\<In the case of vector C cycles reaching 100 or in case you \
want to delete some point(s) in vector C due to a larger than 33% drop/rise \
in vector h, please say in which point you want to stop (number) or, if all \
is ok, write the number 0.\>"]; \ 
    If[lma \[NotEqual] 0, xcycl = {}; l = 1; 
      While[l < lma + 1, xcycl = Append[xcycl, c[\([l]\)]]; l = l + 1]; 
      Print["\<New Vector C:   \>", c = xcycl]]\)], "Input"],

Cell[BoxData[
    \(mma = 
      Input["\<In the case of vector F cycles reaching 100 or in case you \
want to delete some point(s) in vector F due to a larger than 33% drop/rise \
in vector h, please say in which point you want to stop (number) or, if all \
is ok, write the number 0.\>"]; \ 
    If[mma \[NotEqual] 0, ycycl = {}; m = 1; 
      While[m < mma + 1, ycycl = Append[ycycl, f[\([m]\)]]; m = m + 1]; 
      Print["\<New Vector F:   \>", f = ycycl]]\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\( (*\ 
      Routine\ 5\ \((Path)\)\ *) \)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\( (*\ 
      Now\ we\ \ join\ the\ two\ parts\ of\ the\ path; \ 
      c\ includes\ the\ peak, \ 
      f\ starts\ from\ the\ peak\ in\ the\ other\ direction\ but\ does\ not\ \
include\ it; \ 
      we\ first\ create\ g\ that\ is\ f\ reversed\ *) \)\)\)], "Input"],

Cell[BoxData[
    \(compf = Length[f]; k = 0; g = {}; i = 0; 
    While[k < compf, i = f[\([compf - k]\)]; g = Append[g, i]; i = 0; 
      k = k + 1];\)], "Input"],

Cell[BoxData[
    \(\(\( (*\ We\ now\ create\ h, \ 
      that\ combines\ c\ and\ f\ reversed\ \((g)\), \ 
      so\ that\ the\ full\ path\ is\ consistent\ and\ in\ one\ direction\ ; \ 
      it\ is\ important\ to\ show\ vector\ h\ explicitly\ \((for\ eventual\ \
main + lobe\ combinations\ without\ having\ to\ run\ the\ full\ Mathematica\ \
routine\ again; \ Max\ and\ Min\ are\ also\ presented)\); 
      note\ that\ f\ might\ sometimes\ be\ ignored\ \((or\ c)\)\ in\ which\ \
case\ we\ just\ do\ h = 
        c\ \((or\ h = f)\)\ in\ what\ follows\ *) \)\(\ \)\)\)], "Input"],

Cell[BoxData[
    \(\(\(\(h = {}\) \)\(;\)\(\(h = 
        Join[g, c]\) \)\(;\)\(\(Print["\<Vector h:   \>", 
        h]\) \)\(;\)\(\ \)\)\)], "Input"],

Cell[BoxData[
    \(\(\( (*\ now, \ we\ plot\ the\ profile\ of\ the\ path, \ 
      4 - pixel\ per\ 4 - pixel; \ 
      we\ create\ vector\ h3\ which\ will\ be\ very\ relevant\ for\ curve\ \
fitting; \ 
      note\ that\ the\ typical\ profile\ should\ have\ a\ \\\ \(\(or\ /\ 
              or\)\(\ \)\(/\)\)\\\ shape, \ \(\(never\)\(\ \)\(\\/\)\); \ 
      in\ this\ latter\ case\ recreate\ h\ by\ using\ the\ reverse\ order\ in\
\ Join\ \(above : \ Join[c, g]\)\ *) \)\(\[IndentingNewLine]\)\( (*\ 
      In\ addition, \ carefully\ check\ the\ profile, \ 
      since\ it\ must\ follow\ the\ path\ on\ the\ plot; \ 
      when\ the\ path\ is\ too\ vertical, \ this\ is\ easily\ missed; \ so, \ 
      we\ must\ sort\ the\ vector\ \(properly : \ 
          BEFORE\ the\ full\ line\ "\<While...e=Append...l=1;\>"\ \(put : \ "\
\<e=h[[Ordering[h[[All,2]]]]];h=e;e={};\>"\)\)\ *) \)\)\)], "Input"],

Cell[BoxData[
    \( (*\ 
      Now\ is\ also\ \((yet\ another)\)\ moment\ to\ remove\ any\ extremity\ \
points\ that\ drop/rise\ more\ than\ 33  %\ in\ relation\ to\ the\ previous/
          next\ point\ \((use\ routine\ 4\ above\ for\ it, \ 
          as\ we\ have\ done\ in\ this\ example)\)\ *) \)], "Input"],

Cell[BoxData[
    \(l = 1; h3 = {}; n = Length[h]; 
    Print["\<Number of points in vector h:  \>", n]; e = {}; n2 = Length[h]; 
    While[l \[LessEqual] n2, e = Append[e, {h[\([l, 1]\)], h[\([l, 2]\)]}]; 
      l = l + 1]; l = 1; 
    While[l \[LessEqual] n, h3 = Append[h3, h[\([l, 3]\)]]; l = l + 1]; 
    Print["\<Maximum 4-pix of vector h:  \>", Max[h3]]; 
    Print["\<Minimum 4-pix of vector h:  \>", Min[h3]]; 
    ListPlot[h3, PlotRange \[Rule] All, PlotStyle -> PointSize[0.013], 
      Frame -> True, Axes -> False]\)], "Input"],

Cell[BoxData[
    \( (*\ We\ next\ plot\ the\ Pole\ Star\ path\ given\ by\ h {x, 
          y}\ \((x\ is\ each\ first\ coordinate\ h[i, 
              1]\ and\ y\ each\ second\ coordinate\ h[i, 2])\); \ 
      this\ will\ later\ be\ placed\ on\ top\ of\ the\ 3  
        D\ graph\ \((viewed\ as\ a\ 2  D\ contour\ plot; \ 
          presented\ after)\); \ 
      all\ this\ must\ be\ done\ BEFORE\ curve\ fitting\ *) \)], "Input"],

Cell[BoxData[
    \(e1 = 
      ListPlot[e, PlotRange \[Rule] All, 
        PlotStyle -> PointSize[0.013]]\)], "Input"],

Cell[BoxData[
    \(\(\( (*\ 
      First\ the\ 1  x1\ superposition\ *) \)\(\[IndentingNewLine]\)\(e2 = 
      ContourGraphics[a, PlotRange -> {{nx1, nx2}, {ny1, ny2}, {0, topgr}}]; 
    Show[e2]\)\)\)], "Input"],

Cell[BoxData[
    \(Show[e2, e1]\)], "Input"],

Cell[BoxData[
    \( (*\ Now, \ the\ 2  x2\ superposition\ *) \)], "Input"],

Cell[BoxData[
    \(q2 = 
      ContourGraphics[d1, PlotRange -> {{nx1, nx2}, {ny1, ny2}, {0, topgr}}]; 
    Show[q2]\)], "Input"],

Cell[BoxData[
    \(Show[q2, e1]\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\( (*\ 
      Note\ that\ the\ seeing\ determination\ will\ take\ place\ only\ if\ \
there\ is\ no\ "\<LOBE work\>"\ to\ be\ made; \ 
      if\ this\ is\ not\ the\ case, \ 
      jump\ to\ the\ Special\ Lobe\ Section; \ otherwise, \ 
      carry\ on\ with\ Routine\ 6\ *) \)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \( (*\ Now\ begins\ the\ SEEING\ DETERMINATION\ section\ *) \)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\( (*\ 
      Routine\ 6\ \((seeing)\)\ *) \)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \( (*\ Routine\ 6.1\ \((Circle\ Fitting)\)\ *) \)], "Input"],

Cell[BoxData[
    \( (*\ The\ 1  
        st\ attempt\ \((default)\)\ will\ be\ with\ Circle\ Fitting\ *) \)], \
"Input"],

Cell[BoxData[
    \(\[IndentingNewLine]\)], "Input"],

Cell[BoxData[
    \( (*\(\(\ \)\(The\ first\ step\ is\ to\ use\ information\ we\ \(know : \ 
              the\ radius\ of\ the\ circle\); \ 
          this\ is\ the\ exact\ distance\ of\ Polaris\ to\ the\ NCP, \ 
          coming\ from\ the\ exact\ declination\ when\ we\ observed\ Polaris \
. \ Using\ the\ tables\ in\ cadastral . com\ we\ have\ for\ 1  
            st\ Feb\ 2008\ a\ declination\ of\ 89  \.ba18' 
            24.2 "\< and for 5th Mar 2009 89\.ba18'37.2\>" . \ 
              The\ difference\ of\ 13 "\< is irrelevant, giving an error of \
0.5% in the radius, so we pick the mean distance to the NCP (2489.3\>"\)\()\)\
\)\ which, \ 
      divided\ by\ the\ pixel\ size\ of\ 0.77 \*"\"\< gives 3233 pixels \
(rounding to units). \>"*) \)], "Input"],

Cell[BoxData[
    \( (*\ 
      Recalling\ that\ vector\ e\ contains\ the\ data\ in\ the\ format\ we\ \
need, \ \((x, y)\)\ coordinates\ of\ the\ Polaris\ path\ in\ vector\ h, \ 
      we\ then\ fit\ these\ data\ as\ follows\ *) \)], "Input"],

Cell[BoxData[
    \( (*\ It\ is\ important\ to\ modify\ the\ values\ in\ UnitStep, \ 
      since\ these\ focus\ the\ fit\ in\ the\ arc\ section\ where\ the\ data\ \
points\ lie, \ which\ is\ a\ very\ small\ section\ of\ the\ full\ circle; \ 
      we\ use\ the\ \((x1, x2)\), \ \((y1, 
          y2)\)\ limits\ used\ for\ the\ graphs\ above; \ 
      note\ that\ there\ will\ actually\ be\ two\ fits\ \((one\ with\ the\ ' \
+ \(\('\)\(\ \)\(sign\)\); \ 
          the\ other\ with\ the\ ' - \(\('\)\(\ \)\(sign\)\))\); \ 
      sometimes\ one\ fails\ completely\ but, \ more\ typically, \ 
      both\ are\ usually\ similar; \ the\ graphs\ should\ be\ analysed, \ 
      to\ make\ sure\ the\ default\ ' + \(\('\)\(\ \)\(is\)\(\ \)\(always\)\(\
\ \)\(the\)\(\ \)\(best\)\(\ \)\(fit\)\)\ *) \)], "Input"],

Cell[BoxData[
    \( (*\ It\ is\ very\ important\ to\ Sort\ e, \ which\ we\ do\ first; \ 
      note\ that\ if\ the\ path\ is\ too\ vertical, \ 
      Sort\ may\ not\ work\ properly; \ in\ that\ case, \ 
      replace\ "\<Sort[e]\>"\ by\ "\<e[[Ordering[e[[All,2]]]]]\>"\ *) \)], \
"Input"],

Cell[BoxData[
    \(ess = {}; ess = Sort[e]; e = ess;\)], "Input"],

Cell[BoxData[
    \(circfit1 = 
      FindFit[e, 
        atr + Sqrt[
            3233^2 - \((UnitStep[r - nx1]*UnitStep[nx2 - r]*r - 
                    btr)\)^2], {atr, btr}, r, MaxIterations \[Rule] 50000, 
        Method -> LevenbergMarquardt]; atrstr = ToString[circfit1[\([1]\)]]; 
    btrstr = ToString[circfit1[\([2]\)]]; latrstr = StringLength[atrstr]; 
    lbtrstr = StringLength[btrstr]; 
    atr2 = ToExpression[StringTake[atrstr, {8, latrstr}]]; 
    btr2 = ToExpression[StringTake[btrstr, {8, lbtrstr}]];\)], "Input"],

Cell[BoxData[
    \(circfit2 = 
      FindFit[e, 
        atm - Sqrt[
            3233^2 - \((UnitStep[r - nx1]*UnitStep[nx2 - r]*r - 
                    btm)\)^2], {atm, btm}, r, MaxIterations \[Rule] 50000, 
        Method -> LevenbergMarquardt]; atmstr = ToString[circfit2[\([1]\)]]; 
    btmstr = ToString[circfit2[\([2]\)]]; latmstr = StringLength[atmstr]; 
    lbtmstr = StringLength[btmstr]; 
    atm2 = ToExpression[StringTake[atmstr, {8, latmstr}]]; 
    btm2 = ToExpression[StringTake[btmstr, {8, lbtmstr}]];\)], "Input"],

Cell[BoxData[
    \( (*\ We\ now\ build\ each\ graph\ fit\ *) \)], "Input"],

Cell[BoxData[
    \(e10 = 
      Plot[atr2 + Sqrt[3233^2 - \((r - btr2)\)^2], {r, nx1, nx2}, 
        PlotStyle \[Rule] {RGBColor[0, 1, 0], AbsoluteThickness[2]}]; 
    e20 = Plot[atm2 - Sqrt[3233^2 - \((r - btm2)\)^2], {r, nx1, nx2}, 
        PlotStyle \[Rule] {RGBColor[0, 1, 0], AbsoluteThickness[2]}]; 
    Show[e10, e1]; Show[e20, e1]\ ; Show[q2, e10]; Show[q2, e20];\)], "Input"],

Cell[BoxData[
    \( (*\ Finally, \ 
      we\ get\ the\ seeing\ calculations\ for\ the\ circle\ fit; \ 
      we\ pick\ the\ best\ from\ the\ two\ above\ *) \)], "Input"],

Cell[BoxData[
    \(\(ans = 
        Input["\<After careful eye inspection of both plots, in points and in \
greyscale, including curvature analysis, which one is better? If the 1st one \
(or the two cannot be distinguished), please place '1' below (only the \
number). If the 2nd one, please place '2'.\>"];\)\)], "Input"],

Cell[BoxData[
    \(\(If[ans \[Equal] 2, at = atm2; bt = btm2; sinal = \(-1\), at = atr2; 
        bt = btr2; sinal = 1];\)\)], "Input"],

Cell[BoxData[
    \( (*\ 
      We\ use\ the\ distance\ formula\ by\ constructing\ vector\ distancirc\ *) \
\)], "Input"],

Cell[BoxData[
    \(l = 1; n2 = Length[h]; distancirc = {}; somasqr = 0; 
    While[l \[LessEqual] n2, 
      sdf = Minimize[
          Sqrt[\((t - h[\([l, 1]\)])\)^2 + \((at + 
                    sinal*Sqrt[Abs[3233^2 - \((t - bt)\)^2]] - 
                    h[\([l, 2]\)])\)^2], t]; 
      distancirc = Append[distancirc, sdf[\([1]\)]]; 
      somasqr = somasqr + sdf[\([1]\)]^2; l = l + 1]; 
    seeingPh = Median[distancirc]; seeingVis = 2*seeingPh;\)], "Input"],

Cell[BoxData[
    \( (*\ As\ you\ can\ see\ above, \ 
      the\ routine\ just\ calculated\ the\ photographic\ and\ visual\ seeing\ \
values\ from\ the\ mean; \ however, \ 
      we\ will\ print\ out\ \((and\ consider)\)\ only\ the\ photographic\ R . 
          M . S . \ value, \ 
      which\ we\ calculate\ \((and\ print\ out)\)\ next\ \ *) \)], "Input"],

Cell[BoxData[
    \(seeingPh2 = Sqrt[somasqr/n2]; \ seeingVis2 = 2*seeingPh2; 
    sddev = seeingVis2*0.77/
          206265; \ \ Print["\<CCD Seeing (R.M.S.) - CIRCLE FIT in arcseconds \
=  \>", seeingPh2*0.77\ ]\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ Routine\ 6.2\ \((Parabole\ Fitting)\)\ *) \)\)\)], "Input"],

Cell[BoxData[
    \( (*\ For\ comparison, \ 
      we\ have\ actually\ done\ a\ Parabole\ Fit\ too\ \((when\ possible)\), \
\ following\ Moroder\ and\ Righini\ 1973; \ of\ course, \ 
      such\ a\ fit\ is\ the\ only\ one\ available\ when\ Circle\ Fitting\ \
fails\ *) \)], "Input"],

Cell[BoxData[
    \( (*\ The\ parabole\ fit\ is\ simple\ but\ not\ too\ physical; \ and, \ 
      sometimes, \ \(the\ graphs\ are\ just\ stupid\ since\ we\ know\ Polaris\
\ is\ not\ following\ the\ full\ path\ that\ we\ see ... \)\ It\ is, \ 
      thus, \ a\ good\ idea\ to\ always\ get\ the\ alternative\ and\ more\ \
physical\ fit\ of\ a\ circle, \ 
      since\ this\ is\ the\ actual\ motion\ of\ Polaris\ *) \)], "Input"],

Cell[BoxData[{
    \(parabfit1 = FindFit[e, a3 + a2\ xp + a1\ xp^2, \ {a3, a2, a1}, xp]; 
    a3str = ToString[parabfit1[\([1]\)]]; 
    a2str = ToString[parabfit1[\([2]\)]]; 
    a1str = ToString[parabfit1[\([3]\)]]; la3str = StringLength[a3str]; 
    la2str = StringLength[a2str]; la1str = StringLength[a1str]; 
    aa3 = ToExpression[
        StringTake[a3str, {7, la3str}]];\), "\[IndentingNewLine]", 
    \(aa2 = ToExpression[StringTake[a2str, {7, la2str}]]; 
    aa1 = ToExpression[
        StringTake[a1str, {7, la1str}]];\), "\[IndentingNewLine]", 
    \(\(\(e3 = 
      Plot[aa3 + aa2\ xrt + aa1\ xrt^2, {xrt, nx1, nx2}, 
        PlotStyle \[Rule] {RGBColor[0, 1, 0], 
            AbsoluteThickness[2]}]\)\(;\)\(Show[e3, e1]\)\(;\)\(Show[q2, 
      e3]\)\(;\)\(l = 1\)\(;\)\(n2 = 
      Length[e]\)\(;\)\(distancia = {}\)\(;\)\(somasqr = 0\)\(;\)\(While[
      l \[LessEqual] n2, 
      sdf = Minimize[
          Sqrt[\((tpy - e[\([l, 1]\)])\)^2 + \((aa3 + aa2\ tpy + aa1\ tpy^2 - 
                    e[\([l, 2]\)])\)^2], tpy]; 
      distancia = Append[distancia, sdf[\([1]\)]]; 
      somasqr = somasqr + sdf[\([1]\)]^2; l = l + 1]\)\(;\)\(seeingPh = 
      Median[distancia]\)\(;\)\(seeingVis = 2*seeingPh\)\(;\)\(seeingPh2 = 
      Sqrt[somasqr/n2]\)\(;\)\(\ \)\(seeingVis2 = 2*seeingPh2\)\(;\)\(sddev = 
      seeingVis2*0.77/
          206265\)\(;\)\(\ \)\(Print["\<CCD Seeing (R.M.S.) - PARABOLE FIT in \
arcseconds =  \>", 
      seeingPh2*0.77\ ]\)\(;\)\(\[IndentingNewLine]\)\)\)}], "Input"],

Cell[BoxData[
    \( (*\ 
      The\ seeing\ just\ calculated\ is\ the\ transversal\ one; \ \(Moroder\  \
&\)\ Righini\ 1973\ use\ a\ different\ method\ to\ deduce\ the\ total\ \
seeing; \ however, \ 
      they\ used\ photography\ \((so, \ 
          only\ statistically\ could\ they\ evaluate\ the\ deviations\ from\ \
the\ parabole)\), \ we\ use\ a\ CCD, \ 
      much\ more\ precise\ and\ we\ could\ calculate\ the\ deviations\ \
exactly; \ also, \ 
      we\ will\ not\ extrapolate\ to\ large\ telescopes\ as\ they\ did\ \
\((our\ value\ is\ conservative\ anyway, \ 
          this\ way\  - \ 
              it\ would\ only\ get\ even\ smaller ... )\)\ *) \)], "Input"],

Cell[BoxData[
    \( (*\ Still, \ 
      we\ can\ think\ a\ bit\ on\ the\ longitudinal\ seeing\ component; \ 
      this\ is, \ no\ doubt, \ 
      embedded\ in\ the\ brighthness\ variations\ along\ the\ track; \ 
      we\ are\ undersampling\ the\ pole\ star\ "\<Airy Disk\>"\ so\ we\ will\ \
have\ this\ in\ \(mind : \ 
          some\ \((unknown)\)\ part\ of\ the\ brightness\ variations\ is\ \
likely\ due\ to\ intrapixel\ variations\); \ however, \ we\ would\ hope, \ 
      in\ general, \ 
      some\ relation\ between\ bad\ seeing\ and\ brightness\ variations\ \
\((which\ would\ increase)\); \ for\ all\ these, \ at\ this\ stage, \ 
      we\ will\ do\ these\ calculations\ and\ build\ a\ new\ vector\ \
flutuacoes\ with\ all\ brighness\ variations, \ point\ after\ point; \ 
      note\ that\ there\ is\ some\ "\<atenuation\>"\ factor\ to\ bear\ in\ \
mind\ since\ horizontal\ or\ vertical\ motions\ along\ the\ chip\ will\ \
include\ more\ of\ the\ Airy\ disk\ in\ the\ 4 - 
        pixel\ set\ than\ 45\ degree\ movements\ *) \)], "Input"],

Cell[BoxData[
    \( (*\ 
      We\ then\ start\ to\ determine\ the\ angle\ between\ the\ first\ and\ \
last\ point\ of\ h\ and\ the\ consequent\ reduction\ factor\ called\ \
fatorredutor; 
      we\ will\ multiply\ this\ value\ by\ the\ final\ percentage\ \
calculated; \ 
      we\ will\ compare\ each\ of\ the\ ADU\ values\ in\ the\ h\ vector\ \
extremes\ \((the\ lowest\ ones)\)\ with\ the\ peak; \ 
      the\ compensation\ factor\ will\ allow\ datasets\ to\ be\ comparable\ \
in\ the\ overall\ changes\ \((percentages)\); \ 
      this\ is\ also\ an\ issue\ for\ point - to - 
        point\ fluctuations\ as\ we\ go\ along\ the\ vector\ h, \ 
      although\ in\ this\ case\ only\ one\ angle\ out\ of\ horizontal/
          vertical\ is\ \(possible : \ 45\ degrees\)\ *) \)], "Input"],

Cell[BoxData[
    \(u1 = h[\([1, 1]\)]; t1 = h[\([1, 2]\)]; w1 = h[\([n2, 1]\)]; 
    v1 = h[\([n2, 2]\)]; angyy = Abs[ArcTan[Abs[u1 - w1]/Abs[t1 - v1]]]; 
    If[angyy >  .78539816344, angyy = Pi/2 - angyy]; 
    fatorredutor = Sqrt[1 + \((Tan[angyy])\)^2];\)], "Input"],

Cell[BoxData[
    \( (*\ Now\ we\ calculate\ the\ overall\ fluctuation, \ 
      calibrated\ by\ the\ fatorredutor; \ 
      if\ horizontal\ or\ vertical\ fatorredutor = 1; \ 
      note\ that\ the\ typical\ path\ described\ by\ the\ vector\ h\ is\ from\
\ a\ low\ value\ to\ a\ maximum\ value\ to\ a\ low\ value\ again; \ 
      even\ it\ this\ is\ not\ the\ case, \ 
      we\ always\ use\ the\ extremes\ for\ the\ calculation; \ 
      recall\ that\ the\ peak\ value\ \((and\ its\ coordinates)\)\ are\ the\ \
first\ element\ of\ vector\ c\ *) \)], "Input"],

Cell[BoxData[
    \(fluct1 = \((c[\([1, 3]\)] - h[\([1, 3]\)])\)*
        fatorredutor*100/c[\([1, 3]\)]; 
    fluct2 = \((c[\([1, 3]\)] - h[\([n2, 3]\)])\)*
        fatorredutor*100/c[\([1, 3]\)]; \ 
    Print["\<Peak-to-1st-path-point fluctuation: \>", fluct1, \ "\<%\>"]; 
    Print["\<Peak-to-last-path-point fluctuation: \>", 
      fluct2, \ "\<%\>"];\)], "Input"],

Cell[BoxData[
    \( (*\ Now, \ point - to - point\ fluctuations\ accross\ h; \ 
      we\ build\ vector\ flutuacoes; \ 
      for\ simplicity\ we\ use\ the\ same\ angle\ formula\ as\ above, \ 
      otherwise\ we\ would\ have\ to\ have\ a\ set\ of\ IFs; \ 
      the\ final\ value\ will\ be\ a\ median\ of\ all\ values, \ 
      just\ like\ we\ did\ for\ the\ seeing\ \((transversal)\)\ *) \)], \
"Input"],

Cell[BoxData[
    \(l = 1; n2 = Length[h]; flutuacoes = {}; fluct = 0; 
    While[l < n2, u1 = h[\([l, 1]\)]; t1 = h[\([l, 2]\)]; 
      w1 = h[\([l + 1, 1]\)]; v1 = h[\([l + 1, 2]\)]; 
      If[t1 - v1 \[NotEqual] 0, 
        angyy = Abs[ArcTan[Abs[u1 - w1]/Abs[t1 - v1]]]; 
        If[angyy >  .78539816344, angyy = Pi/2 - angyy]; 
        fatorredutor = Sqrt[1 + \((Tan[angyy])\)^2], fatorredutor = 1]; 
      fluct = Abs[\((h[\([l, 3]\)] - h[\([l + 1, 3]\)])\)*fatorredutor]*100/
            h[\([l, 3]\)]; flutuacoes = Append[flutuacoes, fluct]; 
      l = l + 1]; seeinglong = Median[flutuacoes]; 
    Print["\<Longitudinal 'seeing' (percentage of brightness variations along \
path):  \>", seeinglong, \ "\< %\>"]\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ Routine\ 6.3\ \(("\<Horizontal\>"\ Fitting)\)\ *) \)\)\)], "Input"],

Cell[BoxData[
    \( (*\ The\ 3  
        rd\ attempt\ , \ \(when\ Circle\  &\)\ Parabole\ fitting\ both\ fail, \
\ is\ the\ "\<Horizontal Method\>"\ *) \)], "Input"],

Cell[BoxData[
    \(\(\( (*\ 
      In\ the\ rare\ situation\ that\ the\ values\ of\ seeing\ are\ ridiculus\
\ \((usually\ because\ the\ fitted\ arc\ of\ circle\ is\ too\ vertical)\)\ \
then\ run\ the\ following\ routine\ that\ calculates\ the\ horizontal\ \
distance\ of\ each\ point\ of\ vector\ h\ to\ the\ arc\ of\ circle\ *) \)\(\ \
\)\)\)], "Input"],

Cell[BoxData[
    \(\(\(l = 1\) \); \(\(n2 = 
        Length[h]\) \); \(\(distancirc = {}\) \); \(\(xarcocirculo = {}\) \); \
\(\(valor = 0\) \); \(\(sdf = 0\) \); \(\(somasqr = 0\) \); \(\(ArcoCirculo = 
        at + sinal*Sqrt[3233^2 - \((t - bt)\)^2]\) \); \(\(While[
        l \[LessEqual] n2, 
        xarcocirculo = Solve[ArcoCirculo == h[\([l, 2]\)], t]; 
        valor = t\  /. \ xarcocirculo[\([1]\)]; 
        sdf = Abs[valor - h[\([l, 1]\)]]; 
        distancirc = Append[distancirc, sdf]; somasqr = somasqr + sdf^2; 
        l = l + 1; xarcocirculo = {}; sdf = 0; 
        valor = 0]\) \); \(\(seeingPh = 
        Median[distancirc]\) \); \(\(seeingVis = 
        2*seeingPh\) \); \(\(seeingPh2 = 
        Sqrt[somasqr/n2]\) \); \ \(\(seeingVis2 = 
        2*seeingPh2\) \); \(\(sddev = 
        seeingVis2*0.77/
            206265\) \); \ \ Print["\<CCD Seeing (R.M.S.) - HORIZONTAL FIT in \
arcseconds =  \>", seeingPh2*0.77\ ];\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\( (*\ 
      Section\  - \ Lobe\ Fit\ \((if\ necessary)\); \ 
      Lobe\ 1\ *) \)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \( (*\(\(\(\(\ \)\(What\ follows\ essentially\ repeats\ all\ the\ \
routines\ above\ except\ for\ two\ crucial\ \(steps : \ 
                  i\)\)\()\)\)\ \ it\ saves\ the\ current\ vector\ h\ as\ \
vector\ h1\ for\ later\ combination; \ 
          ii\)\()\)\)\ the\ initial\ file\ imported\ \((into\ vector\ a)\)\ \
will\ be\ transformed\ upon\ input\ into\ vector\ al1, \ 
      before\ proceeding\ as\ before; \ 
      some\ common\ values\ to\ the\ above\ routines\ are\ not\ input\ here\ *) \
\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ Routine\ 1\  - \ 
        Initial\ look\ at\ the\ data\ *) \)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[{
    \(h1 = h; \ b = Transpose[a]; 
    ax1 = Input["\<Please insert the limits of the rectangle you want to \
'remove' from the data as follows:  x1?\>"]; ax2 = Input["\<x2?\>"]; 
    ay1 = Input["\<y1?\>"]; 
    ay2 = Input["\<y2?\>"]; \ \ Print["\<You defined the following EXTRACTION \
square: x1= \>", ax1, "\<, x2= \>", ax2, "\<, y1= \>", \ 
      ay1, \ "\<, y2= \>", \ ay2]; x = ax1; y = ay1; 
    While[y < ay2 + 1, \ 
      While[x < ax2 + 1, \ b[\([x, y]\)] = 100; x = x + 1]; x = ax1; 
      y = y + 1];\), "\[IndentingNewLine]", 
    \(al1 = b; m = Max[al1]; 
    Print["\<Peak value of Lobe 1 (ADU):  \>", m];\)}], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ \ Routine\ 2\ is\ skipped\ here\ *) \)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ Routine\ 3\ \((Main)\)\  - \ 
        Organizing\ the\ data\ into\ 2  x2\ pixels, \ 
      preparing\ vectors\ c\ and\ f\ for\ analysis; \ be\ patient, \ 
      it\ might\ take\ a\ few\ \(minutes!\)\ *) \)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[{
    \(x = 1; y = 1; mx = 0; my = 0; 
    While[y < maty + 1, \ 
      While[x < matx + 1, \ 
        If[b[\([x, y]\)] == m, mx = x; my = y; x = matx + 1; y = maty + 1, 
          x = x + 1, Print["\<erro\>"]]]; x = 1; y = y + 1]; \ 
    topgr = m + 5000; x = mx - 1; y = my - 1; i = 0; j = 0; u = 100; 
    While[y < my + 2, 
      While[x < mx + 2, 
        If[b[\([x, y]\)] > u, u = b[\([x, y]\)]; i = x; j = y]; x = x + 1]; 
      y = y + 2; x = mx - 1]; 
    If[b[\([mx - 1, my]\)] > u, u = b[\([mx - 1, my]\)]; i = mx - 1; 
      j = my]; \ 
    If[b[\([mx + 1, my]\)] > u, u = b[\([mx + 1, my]\)]; i = mx + 1; j = my]; 
    m2 = u; \ m2x = i; m2y = j; esq = 0; drt = 0; cim = 0; bxo = 0; dx = 0; 
    dy = 0; dm = 0; es = 0; d = {}; e = {}; i = 0; j = 0; 
    If[Abs[mx - m2x] \[Equal] 1, i = 1, i = 0]; \ 
    If[Abs[my - m2y] \[Equal] 1, j = 1, j = 0]; 
    If[i + j \[Equal] 2, 
      If[mx > m2x && m2y > my, dx = mx - 0.5; dy = my + 0.5; 
        dm = N[Mean[{m, m2, b[\([mx - 1, my]\)], b[\([mx, my + 1]\)]}]]; 
        es = N[StandardDeviation[{m, m2, b[\([mx - 1, my]\)], 
                b[\([mx, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
        e = Append[e, {dx, dy, es}], 
        If[m2x > mx && m2y > my, dx = mx + 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my]\)], b[\([mx, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], 
          If[m2x > mx && my > m2y, dx = mx + 0.5; dy = my - 0.5; 
            dm = N[Mean[{m, m2, b[\([mx + 1, my]\)], b[\([mx, my - 1]\)]}]]; 
            es = N[StandardDeviation[{m, m2, b[\([mx + 1, my]\)], 
                    b[\([mx, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
            e = Append[e, {dx, dy, es}], 
            If[mx > m2x && my > m2y, dx = mx - 0.5; dy = my - 0.5; 
              dm = N[Mean[{m, m2, b[\([mx - 1, my]\)], 
                      b[\([mx, my - 1]\)]}]]; 
              es = N[StandardDeviation[{m, m2, b[\([mx - 1, my]\)], 
                      b[\([mx, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
              e = Append[e, {dx, dy, es}]]]]], 
      If[m2x == mx && m2y > my, 
        esq = Max[b[\([mx - 1, my + 1]\)], b[\([mx - 1, my]\)]]; 
        dir = Max[b[\([mx + 1, my + 1]\)], b[\([mx + 1, my]\)]]; 
        If[esq > dir, dx = mx - 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx - 1, my]\)], 
                  b[\([mx - 1, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx - 1, my]\)], 
                  b[\([mx - 1, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], dx = mx + 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx + 1, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx + 1, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}]]]; 
      If[m2x == mx && m2y < my, 
        esq = Max[b[\([mx - 1, my - 1]\)], b[\([mx - 1, my]\)]]; 
        dir = Max[b[\([mx + 1, my - 1]\)], b[\([mx + 1, my]\)]]; 
        If[esq > dir, dx = mx - 0.5; dy = my - 0.5; 
          dm = N[Mean[{m, m2, b[\([mx - 1, my]\)], 
                  b[\([mx - 1, my - 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx - 1, my]\)], 
                  b[\([mx - 1, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], dx = mx + 0.5; dy = my - 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx + 1, my - 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx + 1, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}]]]; 
      If[m2y == my && m2x < mx, 
        cim = Max[b[\([mx - 1, my + 1]\)], b[\([mx, my + 1]\)]]; 
        bxo = Max[b[\([mx - 1, my - 1]\)], b[\([mx, my - 1]\)]]; 
        If[cim > bxo, dx = mx - 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx - 1, my + 1]\)], 
                  b[\([mx, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx - 1, my + 1]\)], 
                  b[\([mx, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], dx = mx - 0.5; dy = my - 0.5; 
          dm = N[Mean[{m, m2, b[\([mx - 1, my - 1]\)], 
                  b[\([mx, my - 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx - 1, my - 1]\)], 
                  b[\([mx, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}]]]; 
      If[m2y == my && m2x > mx, 
        cim = Max[b[\([mx + 1, my + 1]\)], b[\([mx, my + 1]\)]]; 
        bxo = Max[b[\([mx + 1, my - 1]\)], b[\([mx, my - 1]\)]]; 
        If[cim > bxo, dx = mx + 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my + 1]\)], 
                  b[\([mx, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my + 1]\)], 
                  b[\([mx, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], dx = mx + 0.5; dy = my - 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my - 1]\)], 
                  b[\([mx, my - 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my - 1]\)], 
                  b[\([mx, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}]]]]; mx = dx; my = dy; m = dm; ordd = 0; 
    ul = 0; ur = 0; dl = 0; dr = 0; i = dx - 0.5; j = dy - 0.5; dx = 0; 
    dy = 0; dm = 0; es = 0; xlo = 1; xhi = matx; ylo = 1; yhi = maty; d = {}; 
    e = {}; d1 = {}; \ d2 = {}; \ 
    If[EvenQ[IntegerPart[j]] \[Equal] True, ylo = 2; yhi = maty - 1]; 
    If[EvenQ[IntegerPart[i]] \[Equal] False, xhi = matx - 1, xlo = 2]; 
    y = ylo; v = yhi; x = xlo; l = xhi; 
    While[y < v + 1, 
      While[x < l + 1, dl = b[\([x, y]\)]; dr = b[\([x + 1, y]\)]; 
        ul = b[\([x, y + 1]\)]; ur = b[\([x + 1, y + 1]\)]; dx = x + 0.5; 
        dy = y + 0.5; dm = N[Mean[{dl, dr, ul, ur}]]; 
        es = N[StandardDeviation[{dl, dr, ul, ur}]]; 
        d = Append[d, {dx, dy, dm}]; d2 = Append[d2, dm]; 
        d2 = Append[d2, dm]; 
        If[dx == mx && dy \[Equal] my, ordd = Length[d]]; 
        e = Append[e, {dx, dy, es}]; x = x + 2]; x = xlo; 
      d1 = Append[d1, d2]; d1 = Append[d1, d2]; d2 = {}; y = y + 2]; 
    d3 = Transpose[d1]; c = {{mx, my, m}}; k = 0; l = 0; m2xyADU = 0; 
    u = 10001; ordem = ordd - IntegerPart[\((matx + 3)\)/2]; 
    While[l < 2, 
      While[k < 3, ordem = ordem + 1; 
        If[d[\([ordem, 3]\)] > 10000 && d[\([ordem, 3]\)] > u, 
          u = d[\([ordem, 3]\)]; m2xyADU = ordem]; k = k + 1]; k = 0; 
      ordem = ordem + matx - 4; l = l + 1]; ordem = ordd - 3; 
    While[k < 2, ordem = ordem + 2; 
      If[d[\([ordem, 3]\)] > 10000 && d[\([ordem, 3]\)] > u, 
        u = d[\([ordem, 3]\)]; m2xyADU = ordem]; k = k + 1]; m2 = u; 
    m2x = d[\([m2xyADU, 1]\)]; m2y = d[\([m2xyADU, 2]\)]; 
    If[m2 > 10000, c = Append[c, {m2x, m2y, m2}]]; k = 0; l = 0; m3xyADU = 0; 
    u = 10001; ordem = ordd - IntegerPart[\((matx + 3)\)/2]; f = {}; 
    While[l < 2, 
      While[k < 3, ordem = ordem + 1; 
        If[d[\([ordem, 3]\)] > 10000 && d[\([ordem, 3]\)] > u && 
            d[\([ordem, 1]\)] \[NotEqual] m2x && 
            d[\([ordem, 2]\)] \[NotEqual] m2y, u = d[\([ordem, 3]\)]; 
          m3xyADU = ordem]; k = k + 1]; k = 0; ordem = ordem + matx - 4; 
      l = l + 1]; ordem = ordd - 3; 
    While[k < 2, ordem = ordem + 2; 
      If[d[\([ordem, 3]\)] > 10000 && d[\([ordem, 3]\)] > u && 
          d[\([ordem, 1]\)] \[NotEqual] m2x && 
          d[\([ordem, 2]\)] \[NotEqual] m2y, u = d[\([ordem, 3]\)]; 
        m3xyADU = ordem]; k = k + 1]; m3 = u; m3x = d[\([m3xyADU, 1]\)]; 
    m3y = d[\([m3xyADU, 2]\)]; If[m3 > 10000, f = Append[f, {m3x, m3y, m3}]]; 
    st = 0; temporario = 0; ciclos1 = 0; m1xyADU = ordd; m2ADU = m2xyADU; 
    val382 = IntegerPart[\((matx - 1)\)/2];\), "\[IndentingNewLine]", 
    \(While[st == 0 && ciclos1 < 100, orddif = m2xyADU - ordd; 
      tempor = m2xyADU; u = 0; chave = Abs[val382 - Abs[orddif]]; ordem = 0; 
      teste = 0; 
      If[d[\([ordd + 2*orddif, 3]\)] > 10000 && 
          d[\([ordd + 2*orddif, 3]\)] > u, u = d[\([ordd + 2*orddif, 3]\)]; 
        ordem = ordd + 2*orddif]; \[IndentingNewLine]If\ [chave == 1, 
        teste = orddif - val382; 
        If[teste \[Equal] 1 || teste \[Equal] \(-matx\) + 2, 
          If[d[\([ordd + 2*orddif - 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif - 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif - 1, 3]\)]; 
            ordem = ordd + 2*orddif - 1]; 
          If[d[\([ordd + orddif + 1, 3]\)] > 10000 && 
              d[\([ordd + orddif + 1, 3]\)] > u, 
            u = d[\([ordd + orddif + 1, 3]\)]; ordem = ordd + orddif + 1], 
          If[d[\([ordd + 2*orddif + 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif + 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif + 1, 3]\)]; 
            ordem = ordd + 2*orddif + 1]; 
          If[d[\([ordd + orddif - 1, 3]\)] > 10000 && 
              d[\([ordd + orddif - 1, 3]\)] > u, 
            u = d[\([ordd + orddif - 1, 3]\)]; ordem = ordd + orddif - 1]], 
        If[chave \[Equal] 0, 
          If[d[\([ordd + 2*orddif + 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif + 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif + 1, 3]\)]; 
            ordem = ordd + 2*orddif + 1]; 
          If[d[\([ordd + 2*orddif - 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif - 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif - 1, 3]\)]; 
            ordem = ordd + 2*orddif - 1], 
          If[d[\([tempor + \((val382 + 1)\)*orddif, 3]\)] > 10000 && 
              d[\([tempor + \((val382 + 1)\)*orddif, 3]\)] > u, 
            u = d[\([tempor + \((val382 + 1)\)*orddif, 3]\)]; 
            ordem = tempor + \((val382 + 1)\)*orddif]; 
          If[d[\([tempor - \((val382 - 1)\)*orddif, 3]\)] > 10000 && 
              d[\([tempor - \((val382 - 1)\)*orddif, 3]\)] > u, 
            u = d[\([tempor - \((val382 - 1)\)*orddif, 3]\)]; 
            ordem = tempor - \((val382 - 1)\)*orddif]]]; m2 = u; 
      m2x = d[\([ordem, 1]\)]; m2y = d[\([ordem, 2]\)]; m2xyADU = ordem; 
      ordd = tempor; ciclos1 = ciclos1 + 1; 
      If[m2 > 10000, c = Append[c, {m2x, m2y, m2}], st = 1]]; st = 0; 
    ciclos2 = 0; orddif = 0; temporario = 0; ordd = m1xyADU; 
    m3ADU = m3xyADU;\), "\[IndentingNewLine]", 
    \(While[st == 0 && ciclos2 < 100, orddif = m3xyADU - ordd; 
      temporario = m3xyADU; u = 0; chave = Abs[val382 - Abs[orddif]]; 
      ordem = 0; teste = 0; 
      If[d[\([ordd + 2*orddif, 3]\)] > 10000 && 
          d[\([ordd + 2*orddif, 3]\)] > u, u = d[\([ordd + 2*orddif, 3]\)]; 
        ordem = ordd + 2*orddif]; \[IndentingNewLine]If\ [chave == 1, 
        teste = orddif - val382; 
        If[teste \[Equal] 1 || teste \[Equal] \(-matx\) + 2, 
          If[d[\([ordd + 2*orddif - 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif - 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif - 1, 3]\)]; 
            ordem = ordd + 2*orddif - 1]; 
          If[d[\([ordd + orddif + 1, 3]\)] > 10000 && 
              d[\([ordd + orddif + 1, 3]\)] > u, 
            u = d[\([ordd + orddif + 1, 3]\)]; ordem = ordd + orddif + 1], 
          If[d[\([ordd + 2*orddif + 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif + 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif + 1, 3]\)]; 
            ordem = ordd + 2*orddif + 1]; 
          If[d[\([ordd + orddif - 1, 3]\)] > 10000 && 
              d[\([ordd + orddif - 1, 3]\)] > u, 
            u = d[\([ordd + orddif - 1, 3]\)]; ordem = ordd + orddif - 1]], 
        If[chave \[Equal] 0, 
          If[d[\([ordd + 2*orddif + 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif + 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif + 1, 3]\)]; 
            ordem = ordd + 2*orddif + 1]; 
          If[d[\([ordd + 2*orddif - 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif - 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif - 1, 3]\)]; 
            ordem = ordd + 2*orddif - 1], 
          If[d[\([temporario + \((val382 + 1)\)*orddif, 3]\)] > 10000 && 
              d[\([temporario + \((val382 + 1)\)*orddif, 3]\)] > u, 
            u = d[\([temporario + \((val382 + 1)\)*orddif, 3]\)]; 
            ordem = temporario + \((val382 + 1)\)*orddif]; 
          If[d[\([temporario - \((val382 - 1)\)*orddif, 3]\)] > 10000 && 
              d[\([temporario - \((val382 - 1)\)*orddif, 3]\)] > u, 
            u = d[\([temporario - \((val382 - 1)\)*orddif, 3]\)]; 
            ordem = temporario - \((val382 - 1)\)*orddif]]]; m3 = u; 
      m3x = d[\([ordem, 1]\)]; m3y = d[\([ordem, 2]\)]; m3xyADU = ordem; 
      ordd = temporario; ciclos2 = ciclos2 + 1; 
      If[m3 > 10000, f = Append[f, {m3x, m3y, m3}], st = 1]]; l = 1; e = {}; 
    n2 = Length[c]; 
    While[l \[LessEqual] n2, e = Append[e, {c[\([l, 1]\)], c[\([l, 2]\)]}]; 
      l = l + 1]; 
    ListPlot[e, PlotRange \[Rule] All, PlotStyle -> PointSize[0.013]]; \ 
    Print["\<Cycles of vector c, which includes the peak:  \>", ciclos1 + 1]; 
    Print["\<Coordinates of the peak (x,y):  \>", 
      d[\([m1xyADU, 1]\)], "\<   \>", d[\([m1xyADU, 2]\)]]; 
    Print["\<Coordinates of the 2nd maximum (x,y):  \>", 
      d[\([m2ADU, 1]\)], "\<   \>", d[\([m2ADU, 2]\)]]; 
    Print["\< Vector C:  \>", c]; l = 1; e = {}; n2 = Length[f]; 
    While[l \[LessEqual] n2, e = Append[e, {f[\([l, 1]\)], f[\([l, 2]\)]}]; 
      l = l + 1]; 
    ListPlot[e, PlotRange \[Rule] All, PlotStyle -> PointSize[0.013]]; 
    Print["\<Cycles of vector f, which does not include the peak:  \>", 
      ciclos2]; 
    Print["\<Coordinates of the 3rd maximum (x,y):  \>", 
      d[\([m3ADU, 1]\)], "\<   \>", d[\([m3ADU, 2]\)]]; 
    Print["\< Vector F:  \>", f];\)}], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ Routine\ 4\ \((Cycles)\)\  - \ 
        Eventually\ redefining\ vector\ c\ and/or\ vector\ f; \ 
      if\ this\ was\ the\ case, \ 
      must\ go\ back\ to\ Routine\ 3\ *) \)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \(lma = 
      Input["\<In the case of vector C cycles reaching 101, please say in \
which point you want to stop (number) or, if all is ok, write the number \
0.\>"]; \ If[lma \[NotEqual] 0, xcycl = {}; l = 1; 
      While[l < lma + 1, xcycl = Append[xcycl, c[\([l]\)]]; l = l + 1]; 
      Print["\<New Vector C:   \>", c = xcycl]]; 
    mma = Input["\<In the case of vector F cycles reaching 100, please say in \
which point you want to stop (number) or, if all is ok, write the number \
0.\>"]; \ If[mma \[NotEqual] 0, ycycl = {}; m = 1; 
      While[m < mma + 1, ycycl = Append[ycycl, f[\([m]\)]]; m = m + 1]; 
      Print["\<New Vector F:   \>", f = ycycl]];\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ Routine\ 5\ \((Path)\)\  - \ 
        edit\ "\<h2=Join[g,c]\>"\ if\ necessary; \ and/or, \ if\ necessary, \ 
      BEFORE\ the\ full\ line\ "\<While...e=Append...l=1;\>"\ \(put : \ \
"\<e=h2[[Ordering[h2[[All,2]]]]];h2=e;e={};\>"\)\ *) \)\(\ \)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], \
"Input"],

Cell[BoxData[
    \(compf = Length[f]; k = 0; g = {}; i = 0; 
    While[k < compf, i = f[\([compf - k]\)]; g = Append[g, i]; i = 0; 
      k = k + 1]; h2 = Join[g, c]; Print["\<Vector h:   \>", h2]; l = 1; 
    hl3 = {}; n = Length[h2]; e = {}; n2 = Length[h2]; 
    While[l \[LessEqual] n2, e = Append[e, {h2[\([l, 1]\)], h2[\([l, 2]\)]}]; 
      l = l + 1]; l = 1; 
    While[l \[LessEqual] n, hl3 = Append[hl3, h2[\([l, 3]\)]]; l = l + 1]; 
    ListPlot[hl3, PlotRange \[Rule] All, PlotStyle -> PointSize[0.013], 
      Frame -> True, Axes -> False]; 
    el1 = ListPlot[e, PlotRange \[Rule] All, PlotStyle -> PointSize[0.013]]; 
    el2 = ContourGraphics[a, 
        PlotRange -> {{nx1, nx2}, {ny1, ny2}, {0, topgr}}]; Show[el2]; 
    Show[el2, el1]; 
    ql2 = ContourGraphics[d1, 
        PlotRange -> {{nx1, nx2}, {ny1, ny2}, {0, topgr}}]; Show[ql2]; 
    Show[ql2, el1];\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\(\(Section\  - \ Lobe\ Fit\ \((if\ necessary)\); \ 
    Lobe\ 2\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)
    \)\)\)], "Input"],

Cell[BoxData[
    \( (*\ If\ there\ is\ a\ second\ Lobe, \ this\ section\ must\ be\ run; \ 
      otherwise, \ 
      jump\ to\ the\ "\<Section - combining lobes\>"\ at\ the\ end\ *) \)], \
"Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ Routine\ 1\  - \ 
        Initial\ look\ at\ the\ data\ *) \)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[{
    \(b = Transpose[a]; 
    ax1 = Input["\<Please insert the limits of the rectangle you want to \
'remove' from the data as follows:  x1?\>"]; ax2 = Input["\<x2?\>"]; 
    ay1 = Input["\<y1?\>"]; 
    ay2 = Input["\<y2?\>"]; \ \ Print["\<You defined the following EXTRACTION \
square: x1= \>", ax1, "\<, x2= \>", ax2, "\<, y1= \>", \ 
      ay1, \ "\<, y2= \>", \ ay2]; x = ax1; y = ay1; 
    While[y < ay2 + 1, \ 
      While[x < ax2 + 1, \ b[\([x, y]\)] = 100; x = x + 1]; x = ax1; 
      y = y + 1];\), "\[IndentingNewLine]", 
    \(al2 = b; m = Max[al2]; 
    Print["\<Peak value of Lobe 2 (ADU):  \>", m];\)}], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\( (*\ 
      Routine\ 2\ is\ skipped\ here\ *) \)\)\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ Routine\ 3\ \((Main)\)\  - \ 
        Organizing\ the\ data\ into\ 2  x2\ pixels, \ 
      preparing\ vectors\ c\ and\ f\ for\ analysis; \ be\ patient, \ 
      it\ might\ take\ a\ few\ \(minutes!\)\ *) \)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[{
    \(x = 1; y = 1; mx = 0; my = 0; 
    While[y < maty + 1, \ 
      While[x < matx + 1, \ 
        If[b[\([x, y]\)] == m, mx = x; my = y; x = matx + 1; y = maty + 1, 
          x = x + 1, Print["\<erro\>"]]]; x = 1; y = y + 1]; \ 
    topgr = m + 5000; x = mx - 1; y = my - 1; i = 0; j = 0; u = 100; 
    While[y < my + 2, 
      While[x < mx + 2, 
        If[b[\([x, y]\)] > u, u = b[\([x, y]\)]; i = x; j = y]; x = x + 1]; 
      y = y + 2; x = mx - 1]; 
    If[b[\([mx - 1, my]\)] > u, u = b[\([mx - 1, my]\)]; i = mx - 1; 
      j = my]; \ 
    If[b[\([mx + 1, my]\)] > u, u = b[\([mx + 1, my]\)]; i = mx + 1; j = my]; 
    m2 = u; \ m2x = i; m2y = j; esq = 0; drt = 0; cim = 0; bxo = 0; dx = 0; 
    dy = 0; dm = 0; es = 0; d = {}; e = {}; i = 0; j = 0; 
    If[Abs[mx - m2x] \[Equal] 1, i = 1, i = 0]; \ 
    If[Abs[my - m2y] \[Equal] 1, j = 1, j = 0]; 
    If[i + j \[Equal] 2, 
      If[mx > m2x && m2y > my, dx = mx - 0.5; dy = my + 0.5; 
        dm = N[Mean[{m, m2, b[\([mx - 1, my]\)], b[\([mx, my + 1]\)]}]]; 
        es = N[StandardDeviation[{m, m2, b[\([mx - 1, my]\)], 
                b[\([mx, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
        e = Append[e, {dx, dy, es}], 
        If[m2x > mx && m2y > my, dx = mx + 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my]\)], b[\([mx, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], 
          If[m2x > mx && my > m2y, dx = mx + 0.5; dy = my - 0.5; 
            dm = N[Mean[{m, m2, b[\([mx + 1, my]\)], b[\([mx, my - 1]\)]}]]; 
            es = N[StandardDeviation[{m, m2, b[\([mx + 1, my]\)], 
                    b[\([mx, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
            e = Append[e, {dx, dy, es}], 
            If[mx > m2x && my > m2y, dx = mx - 0.5; dy = my - 0.5; 
              dm = N[Mean[{m, m2, b[\([mx - 1, my]\)], 
                      b[\([mx, my - 1]\)]}]]; 
              es = N[StandardDeviation[{m, m2, b[\([mx - 1, my]\)], 
                      b[\([mx, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
              e = Append[e, {dx, dy, es}]]]]], 
      If[m2x == mx && m2y > my, 
        esq = Max[b[\([mx - 1, my + 1]\)], b[\([mx - 1, my]\)]]; 
        dir = Max[b[\([mx + 1, my + 1]\)], b[\([mx + 1, my]\)]]; 
        If[esq > dir, dx = mx - 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx - 1, my]\)], 
                  b[\([mx - 1, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx - 1, my]\)], 
                  b[\([mx - 1, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], dx = mx + 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx + 1, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx + 1, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}]]]; 
      If[m2x == mx && m2y < my, 
        esq = Max[b[\([mx - 1, my - 1]\)], b[\([mx - 1, my]\)]]; 
        dir = Max[b[\([mx + 1, my - 1]\)], b[\([mx + 1, my]\)]]; 
        If[esq > dir, dx = mx - 0.5; dy = my - 0.5; 
          dm = N[Mean[{m, m2, b[\([mx - 1, my]\)], 
                  b[\([mx - 1, my - 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx - 1, my]\)], 
                  b[\([mx - 1, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], dx = mx + 0.5; dy = my - 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx + 1, my - 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my]\)], 
                  b[\([mx + 1, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}]]]; 
      If[m2y == my && m2x < mx, 
        cim = Max[b[\([mx - 1, my + 1]\)], b[\([mx, my + 1]\)]]; 
        bxo = Max[b[\([mx - 1, my - 1]\)], b[\([mx, my - 1]\)]]; 
        If[cim > bxo, dx = mx - 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx - 1, my + 1]\)], 
                  b[\([mx, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx - 1, my + 1]\)], 
                  b[\([mx, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], dx = mx - 0.5; dy = my - 0.5; 
          dm = N[Mean[{m, m2, b[\([mx - 1, my - 1]\)], 
                  b[\([mx, my - 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx - 1, my - 1]\)], 
                  b[\([mx, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}]]]; 
      If[m2y == my && m2x > mx, 
        cim = Max[b[\([mx + 1, my + 1]\)], b[\([mx, my + 1]\)]]; 
        bxo = Max[b[\([mx + 1, my - 1]\)], b[\([mx, my - 1]\)]]; 
        If[cim > bxo, dx = mx + 0.5; dy = my + 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my + 1]\)], 
                  b[\([mx, my + 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my + 1]\)], 
                  b[\([mx, my + 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}], dx = mx + 0.5; dy = my - 0.5; 
          dm = N[Mean[{m, m2, b[\([mx + 1, my - 1]\)], 
                  b[\([mx, my - 1]\)]}]]; 
          es = N[StandardDeviation[{m, m2, b[\([mx + 1, my - 1]\)], 
                  b[\([mx, my - 1]\)]}]]; d = Append[d, {dx, dy, dm}]; 
          e = Append[e, {dx, dy, es}]]]]; mx = dx; my = dy; m = dm; ordd = 0; 
    ul = 0; ur = 0; dl = 0; dr = 0; i = dx - 0.5; j = dy - 0.5; dx = 0; 
    dy = 0; dm = 0; es = 0; xlo = 1; xhi = matx; ylo = 1; yhi = maty; d = {}; 
    e = {}; d1 = {}; \ d2 = {}; \ 
    If[EvenQ[IntegerPart[j]] \[Equal] True, ylo = 2; yhi = maty - 1]; 
    If[EvenQ[IntegerPart[i]] \[Equal] False, xhi = matx - 1, xlo = 2]; 
    y = ylo; v = yhi; x = xlo; l = xhi; 
    While[y < v + 1, 
      While[x < l + 1, dl = b[\([x, y]\)]; dr = b[\([x + 1, y]\)]; 
        ul = b[\([x, y + 1]\)]; ur = b[\([x + 1, y + 1]\)]; dx = x + 0.5; 
        dy = y + 0.5; dm = N[Mean[{dl, dr, ul, ur}]]; 
        es = N[StandardDeviation[{dl, dr, ul, ur}]]; 
        d = Append[d, {dx, dy, dm}]; d2 = Append[d2, dm]; 
        d2 = Append[d2, dm]; 
        If[dx == mx && dy \[Equal] my, ordd = Length[d]]; 
        e = Append[e, {dx, dy, es}]; x = x + 2]; x = xlo; 
      d1 = Append[d1, d2]; d1 = Append[d1, d2]; d2 = {}; y = y + 2]; 
    d3 = Transpose[d1]; c = {{mx, my, m}}; k = 0; l = 0; m2xyADU = 0; 
    u = 10001; ordem = ordd - IntegerPart[\((matx + 3)\)/2]; 
    While[l < 2, 
      While[k < 3, ordem = ordem + 1; 
        If[d[\([ordem, 3]\)] > 10000 && d[\([ordem, 3]\)] > u, 
          u = d[\([ordem, 3]\)]; m2xyADU = ordem]; k = k + 1]; k = 0; 
      ordem = ordem + matx - 4; l = l + 1]; ordem = ordd - 3; 
    While[k < 2, ordem = ordem + 2; 
      If[d[\([ordem, 3]\)] > 10000 && d[\([ordem, 3]\)] > u, 
        u = d[\([ordem, 3]\)]; m2xyADU = ordem]; k = k + 1]; m2 = u; 
    m2x = d[\([m2xyADU, 1]\)]; m2y = d[\([m2xyADU, 2]\)]; 
    If[m2 > 10000, c = Append[c, {m2x, m2y, m2}]]; k = 0; l = 0; m3xyADU = 0; 
    u = 10001; ordem = ordd - IntegerPart[\((matx + 3)\)/2]; f = {}; 
    While[l < 2, 
      While[k < 3, ordem = ordem + 1; 
        If[d[\([ordem, 3]\)] > 10000 && d[\([ordem, 3]\)] > u && 
            d[\([ordem, 1]\)] \[NotEqual] m2x && 
            d[\([ordem, 2]\)] \[NotEqual] m2y, u = d[\([ordem, 3]\)]; 
          m3xyADU = ordem]; k = k + 1]; k = 0; ordem = ordem + matx - 4; 
      l = l + 1]; ordem = ordd - 3; 
    While[k < 2, ordem = ordem + 2; 
      If[d[\([ordem, 3]\)] > 10000 && d[\([ordem, 3]\)] > u && 
          d[\([ordem, 1]\)] \[NotEqual] m2x && 
          d[\([ordem, 2]\)] \[NotEqual] m2y, u = d[\([ordem, 3]\)]; 
        m3xyADU = ordem]; k = k + 1]; m3 = u; m3x = d[\([m3xyADU, 1]\)]; 
    m3y = d[\([m3xyADU, 2]\)]; If[m3 > 10000, f = Append[f, {m3x, m3y, m3}]]; 
    st = 0; temporario = 0; ciclos1 = 0; m1xyADU = ordd; m2ADU = m2xyADU; 
    val382 = IntegerPart[\((matx - 1)\)/2];\), "\[IndentingNewLine]", 
    \(While[st == 0 && ciclos1 < 100, orddif = m2xyADU - ordd; 
      tempor = m2xyADU; u = 0; chave = Abs[val382 - Abs[orddif]]; ordem = 0; 
      teste = 0; 
      If[d[\([ordd + 2*orddif, 3]\)] > 10000 && 
          d[\([ordd + 2*orddif, 3]\)] > u, u = d[\([ordd + 2*orddif, 3]\)]; 
        ordem = ordd + 2*orddif]; \[IndentingNewLine]If\ [chave == 1, 
        teste = orddif - val382; 
        If[teste \[Equal] 1 || teste \[Equal] \(-matx\) + 2, 
          If[d[\([ordd + 2*orddif - 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif - 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif - 1, 3]\)]; 
            ordem = ordd + 2*orddif - 1]; 
          If[d[\([ordd + orddif + 1, 3]\)] > 10000 && 
              d[\([ordd + orddif + 1, 3]\)] > u, 
            u = d[\([ordd + orddif + 1, 3]\)]; ordem = ordd + orddif + 1], 
          If[d[\([ordd + 2*orddif + 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif + 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif + 1, 3]\)]; 
            ordem = ordd + 2*orddif + 1]; 
          If[d[\([ordd + orddif - 1, 3]\)] > 10000 && 
              d[\([ordd + orddif - 1, 3]\)] > u, 
            u = d[\([ordd + orddif - 1, 3]\)]; ordem = ordd + orddif - 1]], 
        If[chave \[Equal] 0, 
          If[d[\([ordd + 2*orddif + 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif + 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif + 1, 3]\)]; 
            ordem = ordd + 2*orddif + 1]; 
          If[d[\([ordd + 2*orddif - 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif - 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif - 1, 3]\)]; 
            ordem = ordd + 2*orddif - 1], 
          If[d[\([tempor + \((val382 + 1)\)*orddif, 3]\)] > 10000 && 
              d[\([tempor + \((val382 + 1)\)*orddif, 3]\)] > u, 
            u = d[\([tempor + \((val382 + 1)\)*orddif, 3]\)]; 
            ordem = tempor + \((val382 + 1)\)*orddif]; 
          If[d[\([tempor - \((val382 - 1)\)*orddif, 3]\)] > 10000 && 
              d[\([tempor - \((val382 - 1)\)*orddif, 3]\)] > u, 
            u = d[\([tempor - \((val382 - 1)\)*orddif, 3]\)]; 
            ordem = tempor - \((val382 - 1)\)*orddif]]]; m2 = u; 
      m2x = d[\([ordem, 1]\)]; m2y = d[\([ordem, 2]\)]; m2xyADU = ordem; 
      ordd = tempor; ciclos1 = ciclos1 + 1; 
      If[m2 > 10000, c = Append[c, {m2x, m2y, m2}], st = 1]]; st = 0; 
    ciclos2 = 0; orddif = 0; temporario = 0; ordd = m1xyADU; 
    m3ADU = m3xyADU;\), "\[IndentingNewLine]", 
    \(While[st == 0 && ciclos2 < 100, orddif = m3xyADU - ordd; 
      temporario = m3xyADU; u = 0; chave = Abs[val382 - Abs[orddif]]; 
      ordem = 0; teste = 0; 
      If[d[\([ordd + 2*orddif, 3]\)] > 10000 && 
          d[\([ordd + 2*orddif, 3]\)] > u, u = d[\([ordd + 2*orddif, 3]\)]; 
        ordem = ordd + 2*orddif]; \[IndentingNewLine]If\ [chave == 1, 
        teste = orddif - val382; 
        If[teste \[Equal] 1 || teste \[Equal] \(-matx\) + 2, 
          If[d[\([ordd + 2*orddif - 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif - 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif - 1, 3]\)]; 
            ordem = ordd + 2*orddif - 1]; 
          If[d[\([ordd + orddif + 1, 3]\)] > 10000 && 
              d[\([ordd + orddif + 1, 3]\)] > u, 
            u = d[\([ordd + orddif + 1, 3]\)]; ordem = ordd + orddif + 1], 
          If[d[\([ordd + 2*orddif + 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif + 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif + 1, 3]\)]; 
            ordem = ordd + 2*orddif + 1]; 
          If[d[\([ordd + orddif - 1, 3]\)] > 10000 && 
              d[\([ordd + orddif - 1, 3]\)] > u, 
            u = d[\([ordd + orddif - 1, 3]\)]; ordem = ordd + orddif - 1]], 
        If[chave \[Equal] 0, 
          If[d[\([ordd + 2*orddif + 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif + 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif + 1, 3]\)]; 
            ordem = ordd + 2*orddif + 1]; 
          If[d[\([ordd + 2*orddif - 1, 3]\)] > 10000 && 
              d[\([ordd + 2*orddif - 1, 3]\)] > u, 
            u = d[\([ordd + 2*orddif - 1, 3]\)]; 
            ordem = ordd + 2*orddif - 1], 
          If[d[\([temporario + \((val382 + 1)\)*orddif, 3]\)] > 10000 && 
              d[\([temporario + \((val382 + 1)\)*orddif, 3]\)] > u, 
            u = d[\([temporario + \((val382 + 1)\)*orddif, 3]\)]; 
            ordem = temporario + \((val382 + 1)\)*orddif]; 
          If[d[\([temporario - \((val382 - 1)\)*orddif, 3]\)] > 10000 && 
              d[\([temporario - \((val382 - 1)\)*orddif, 3]\)] > u, 
            u = d[\([temporario - \((val382 - 1)\)*orddif, 3]\)]; 
            ordem = temporario - \((val382 - 1)\)*orddif]]]; m3 = u; 
      m3x = d[\([ordem, 1]\)]; m3y = d[\([ordem, 2]\)]; m3xyADU = ordem; 
      ordd = temporario; ciclos2 = ciclos2 + 1; 
      If[m3 > 10000, f = Append[f, {m3x, m3y, m3}], st = 1]]; l = 1; e = {}; 
    n2 = Length[c]; 
    While[l \[LessEqual] n2, e = Append[e, {c[\([l, 1]\)], c[\([l, 2]\)]}]; 
      l = l + 1]; 
    ListPlot[e, PlotRange \[Rule] All, PlotStyle -> PointSize[0.013]]; \ 
    Print["\<Cycles of vector c, which includes the peak:  \>", ciclos1 + 1]; 
    Print["\<Coordinates of the peak (x,y):  \>", 
      d[\([m1xyADU, 1]\)], "\<   \>", d[\([m1xyADU, 2]\)]]; 
    Print["\<Coordinates of the 2nd maximum (x,y):  \>", 
      d[\([m2ADU, 1]\)], "\<   \>", d[\([m2ADU, 2]\)]]; 
    Print["\< Vector C:  \>", c]; l = 1; e = {}; n2 = Length[f]; 
    While[l \[LessEqual] n2, e = Append[e, {f[\([l, 1]\)], f[\([l, 2]\)]}]; 
      l = l + 1]; 
    ListPlot[e, PlotRange \[Rule] All, PlotStyle -> PointSize[0.013]]; 
    Print["\<Cycles of vector f, which does not include the peak:  \>", 
      ciclos2]; 
    Print["\<Coordinates of the 3rd maximum (x,y):  \>", 
      d[\([m3ADU, 1]\)], "\<   \>", d[\([m3ADU, 2]\)]]; 
    Print["\< Vector F:  \>", f];\)}], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ Routine\ 4\ \((Cycles)\)\  - \ 
        Eventually\ redefining\ vector\ c\ and/or\ vector\ f; \ 
      if\ this\ was\ the\ case, \ 
      must\ go\ back\ to\ Routine\ 3\ *) \)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \(lma = 
      Input["\<In the case of vector C cycles reaching 101, please say in \
which point you want to stop (number) or, if all is ok, write the number \
0.\>"]; \ If[lma \[NotEqual] 0, xcycl = {}; l = 1; 
      While[l < lma + 1, xcycl = Append[xcycl, c[\([l]\)]]; l = l + 1]; 
      Print["\<New Vector C:   \>", c = xcycl]]; 
    mma = Input["\<In the case of vector F cycles reaching 100, please say in \
which point you want to stop (number) or, if all is ok, write the number \
0.\>"]; \ If[mma \[NotEqual] 0, ycycl = {}; m = 1; 
      While[m < mma + 1, ycycl = Append[ycycl, f[\([m]\)]]; m = m + 1]; 
      Print["\<New Vector F:   \>", f = ycycl]];\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ Routine\ 5\ \((Path)\)\  - \ 
        edit\ "\<h4=Join[g,c]\>"\ if\ necessary; \ and/or, \ if\ necessary, \ 
      BEFORE\ the\ full\ line\ "\<While...e=Append...l=1;\>"\ \(put : \ \
"\<e=h4[[Ordering[h4[[All,2]]]]];h4=e;e={};\>"\)\ *) \)\(\ \)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], \
"Input"],

Cell[BoxData[
    \(compf = Length[f]; k = 0; g = {}; i = 0; 
    While[k < compf, i = f[\([compf - k]\)]; g = Append[g, i]; i = 0; 
      k = k + 1]; h4 = Join[g, c]; Print["\<Vector h:   \>", h4]; l = 1; 
    hl3 = {}; n = Length[h4]; e = {}; n2 = Length[h4]; 
    While[l \[LessEqual] n2, e = Append[e, {h4[\([l, 1]\)], h4[\([l, 2]\)]}]; 
      l = l + 1]; l = 1; 
    While[l \[LessEqual] n, hl3 = Append[hl3, h4[\([l, 3]\)]]; l = l + 1]; 
    ListPlot[hl3, PlotRange \[Rule] All, PlotStyle -> PointSize[0.013], 
      Frame -> True, Axes -> False]; 
    el1 = ListPlot[e, PlotRange \[Rule] All, PlotStyle -> PointSize[0.013]]; 
    el2 = ContourGraphics[a, 
        PlotRange -> {{nx1, nx2}, {ny1, ny2}, {0, topgr}}]; Show[el2]; 
    Show[el2, el1]; 
    ql2 = ContourGraphics[d1, 
        PlotRange -> {{nx1, nx2}, {ny1, ny2}, {0, topgr}}]; Show[ql2]; 
    Show[ql2, el1];\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\( (*\ Now, \ 
      depending\ on\ how\ many\ lobes\ \((one\ or\ two)\)\ to\ combine, \ 
      we\ use\ the\ next\ section\ routine\ *) \)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ Section\  - \ 
        combining\ lobes\ *) \)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\(\[IndentingNewLine]\)\)\)], "Input"],

Cell[BoxData[
    \( (*\ We\ combine\ all\ parts\ in\ what\ follows; \ 
      note\ that\ if\ more\ than\ one\ lobe\ then\ "\<hr=Join[h1,h2]\>" 
          must\ be\ changed\ accordingly\ \((e . 
              g . \ "\<hr=Join[h1,h2,h4]\>")\)\  - \ never\ use\ h3; \ 
      if\ more\ than\ two\ lobes\ use\ h5, h6, \ \(\(etc\)\(.\)\); 
      Sort\ is\ a\ crucial\ command\ on\ the\ routine\ so, \ 
      in\ situations\ where\ the\ path\ is\ clearly\ more\ vertical\ than\ \
horizontal, \ Sort\ WILL\ NOT\ work; \ 
      instead\ of\ h = Sort[hr]\ use\ "\<h=hr[[Ordering[hr[[All,2]]]]]\>"; \ 
      finally, \ 
      adjust\ AspectRatio\ to\ make\ a\ realistic\ plot\ of\ the\ path\ *) \
\)], "Input"],

Cell[BoxData[
    \( (*\ Important : \ 
        in\ the\ event\ of\ needing\ to\ change\ vector\ h\ \((remove\ \
points)\)\ either\ because\ of\ a\  > 
          33  %\ rise/drop\ or\ of\ the\ "\<border effect\>", \ 
      you\ will\ have\ to\ do\ so\ by\ hand\ [
          unless\ you\ prefer\ to\ go\ up\ again\ and\ remove\ on\ C\ and\ F\ \
routine \((s)\)], \ 
      simply\ copying\ and\ pasting\ the\ current\ h\ and\ then\ selecting\ \
the\ points\ you\ want\ to\ keep\ to\ define\ a\ new\ vector\ h; \ 
      after\ that, \ simply\ run\ the\ routine\ again, \ 
      AFTER\ removing\ the\ first\ two\ commands\ \((Join\ and\ Sort)\)\ *) \
\)], "Input"],

Cell[BoxData[
    \(hr = Join[h1, h2]; \ h = Sort[hr]; 
    Print["\<The combined and sorted vector H:  \>", h]; l = 1; w3 = {}; 
    e = {}; n = Length[h]; 
    While[l \[LessEqual] n, w3 = Append[w3, h[\([l, 3]\)]]; 
      e = Append[e, {h[\([l, 1]\)], h[\([l, 2]\)]}]; l = l + 1]; 
    ListPlot[w3, PlotRange \[Rule] All, PlotStyle -> PointSize[0.013], 
      Frame \[Rule] True, Axes \[Rule] False]; 
    ec1 = ListPlot[e, PlotRange \[Rule] All, PlotStyle -> PointSize[0.013], 
        AspectRatio \[Rule] 1]; Show[e2, ec1]; Show[q2, ec1]; 
    Print["\<Maximum 4-pix of vector h:  \>", Max[w3]]; 
    Print["\<Minimum 4-pix of vector h:  \>", Min[w3]]; 
    Print["\<Number of points in vector h:  \>", n]; e1 = ec1;\)], "Input"],

Cell[BoxData[
    \(\(\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\
\( (*\ And\ now\ you\ are\ ready\ to\ do\ the\ seeing\ fit\ in\ Routine\ 6\ \
above\ *) \)\(\[IndentingNewLine]\)\(\[IndentingNewLine]\)\(\
\[IndentingNewLine]\)\)\)], "Input"]
},
FrontEndVersion->"5.0 for Microsoft Windows",
ScreenRectangle->{{0, 1920}, {0, 982}},
WindowSize->{674, 766},
WindowMargins->{{Automatic, 434}, {Automatic, 30}}
]

(*******************************************************************
Cached data follows.  If you edit this Notebook file directly, not
using Mathematica, you must remove the line containing CacheID at
the top of  the file.  The cache data will then be recreated when
you save this file from within Mathematica.
*******************************************************************)

(*CellTagsOutline
CellTagsIndex->{}
*)

(*CellTagsIndex
CellTagsIndex->{}
*)

(*NotebookFileOutline
Notebook[{
Cell[1754, 51, 289, 5, 150, "Input"],
Cell[2046, 58, 355, 7, 90, "Input"],
Cell[2404, 67, 434, 8, 90, "Input"],
Cell[2841, 77, 292, 6, 70, "Input"],
Cell[3136, 85, 247, 4, 150, "Input"],
Cell[3386, 91, 557, 10, 130, "Input"],
Cell[3946, 103, 187, 3, 50, "Input"],
Cell[4136, 108, 75, 2, 29, "Input"],
Cell[4214, 112, 209, 4, 50, "Input"],
Cell[4426, 118, 387, 7, 110, "Input"],
Cell[4816, 127, 301, 5, 150, "Input"],
Cell[5120, 134, 536, 8, 150, "Input"],
Cell[5659, 144, 223, 3, 70, "Input"],
Cell[5885, 149, 232, 4, 50, "Input"],
Cell[6120, 155, 163, 3, 30, "Input"],
Cell[6286, 160, 391, 6, 190, "Input"],
Cell[6680, 168, 303, 5, 70, "Input"],
Cell[6986, 175, 303, 6, 110, "Input"],
Cell[7292, 183, 510, 8, 130, "Input"],
Cell[7805, 193, 84, 1, 30, "Input"],
Cell[7892, 196, 435, 9, 110, "Input"],
Cell[8330, 207, 1211, 21, 290, "Input"],
Cell[9544, 230, 4817, 82, 1090, "Input"],
Cell[14364, 314, 841, 12, 250, "Input"],
Cell[15208, 328, 1010, 17, 250, "Input"],
Cell[16221, 347, 281, 5, 110, "Input"],
Cell[16505, 354, 407, 7, 90, "Input"],
Cell[16915, 363, 675, 12, 190, "Input"],
Cell[17593, 377, 860, 16, 210, "Input"],
Cell[18456, 395, 1142, 20, 290, "Input"],
Cell[19601, 417, 2615, 44, 570, "Input"],
Cell[22219, 463, 413, 6, 90, "Input"],
Cell[22635, 471, 2627, 44, 590, "Input"],
Cell[25265, 517, 747, 13, 190, "Input"],
Cell[26015, 532, 494, 9, 130, "Input"],
Cell[26512, 543, 380, 7, 110, "Input"],
Cell[26895, 552, 371, 6, 110, "Input"],
Cell[27269, 560, 332, 5, 90, "Input"],
Cell[27604, 567, 203, 4, 50, "Input"],
Cell[27810, 573, 626, 10, 170, "Input"],
Cell[28439, 585, 83, 1, 70, "Input"],
Cell[28525, 588, 470, 9, 170, "Input"],
Cell[28998, 599, 473, 8, 150, "Input"],
Cell[29474, 609, 473, 8, 170, "Input"],
Cell[29950, 619, 424, 8, 150, "Input"],
Cell[30377, 629, 163, 3, 50, "Input"],
Cell[30543, 634, 578, 9, 150, "Input"],
Cell[31124, 645, 151, 3, 30, "Input"],
Cell[31278, 650, 895, 14, 230, "Input"],
Cell[32176, 666, 312, 5, 70, "Input"],
Cell[32491, 673, 540, 9, 150, "Input"],
Cell[33034, 684, 430, 7, 110, "Input"],
Cell[33467, 693, 119, 3, 30, "Input"],
Cell[33589, 698, 213, 4, 70, "Input"],
Cell[33805, 704, 45, 1, 30, "Input"],
Cell[33853, 707, 75, 1, 30, "Input"],
Cell[33931, 710, 130, 3, 50, "Input"],
Cell[34064, 715, 45, 1, 30, "Input"],
Cell[34112, 718, 131, 2, 110, "Input"],
Cell[34246, 722, 428, 7, 170, "Input"],
Cell[34677, 731, 91, 1, 30, "Input"],
Cell[34771, 734, 128, 2, 70, "Input"],
Cell[34902, 738, 78, 1, 30, "Input"],
Cell[34983, 741, 121, 3, 30, "Input"],
Cell[35107, 746, 52, 1, 50, "Input"],
Cell[35162, 749, 766, 12, 170, "Input"],
Cell[35931, 763, 242, 4, 70, "Input"],
Cell[36176, 769, 803, 12, 190, "Input"],
Cell[36982, 783, 289, 5, 70, "Input"],
Cell[37274, 790, 66, 1, 30, "Input"],
Cell[37343, 793, 532, 10, 150, "Input"],
Cell[37878, 805, 532, 10, 150, "Input"],
Cell[38413, 817, 75, 1, 30, "Input"],
Cell[38491, 820, 385, 6, 130, "Input"],
Cell[38879, 828, 171, 3, 50, "Input"],
Cell[39053, 833, 323, 5, 110, "Input"],
Cell[39379, 840, 136, 2, 30, "Input"],
Cell[39518, 844, 121, 3, 30, "Input"],
Cell[39642, 849, 468, 9, 170, "Input"],
Cell[40113, 860, 357, 6, 110, "Input"],
Cell[40473, 868, 226, 4, 70, "Input"],
Cell[40702, 874, 159, 2, 90, "Input"],
Cell[40864, 878, 282, 5, 70, "Input"],
Cell[41149, 885, 426, 6, 110, "Input"],
Cell[41578, 893, 1511, 28, 370, "Input"],
Cell[43092, 923, 674, 13, 190, "Input"],
Cell[43769, 938, 1051, 17, 250, "Input"],
Cell[44823, 957, 791, 14, 230, "Input"],
Cell[45617, 973, 271, 4, 90, "Input"],
Cell[45891, 979, 559, 9, 150, "Input"],
Cell[46453, 990, 370, 7, 90, "Input"],
Cell[46826, 999, 406, 7, 110, "Input"],
Cell[47235, 1008, 733, 12, 190, "Input"],
Cell[47971, 1022, 167, 2, 90, "Input"],
Cell[48141, 1026, 166, 3, 50, "Input"],
Cell[48310, 1031, 355, 6, 90, "Input"],
Cell[48668, 1039, 955, 18, 210, "Input"],
Cell[49626, 1059, 402, 6, 270, "Input"],
Cell[50031, 1067, 531, 9, 110, "Input"],
Cell[50565, 1078, 247, 4, 150, "Input"],
Cell[50815, 1084, 651, 12, 170, "Input"],
Cell[51469, 1098, 224, 3, 150, "Input"],
Cell[51696, 1103, 391, 6, 190, "Input"],
Cell[52090, 1111, 13918, 243, 3190, "Input"],
Cell[66011, 1356, 361, 6, 170, "Input"],
Cell[66375, 1364, 694, 11, 270, "Input"],
Cell[67072, 1377, 432, 7, 190, "Input"],
Cell[67507, 1386, 891, 16, 250, "Input"],
Cell[68401, 1404, 252, 5, 150, "Input"],
Cell[68656, 1411, 198, 4, 50, "Input"],
Cell[68857, 1417, 247, 4, 150, "Input"],
Cell[69107, 1423, 641, 12, 170, "Input"],
Cell[69751, 1437, 133, 2, 70, "Input"],
Cell[69887, 1441, 391, 6, 190, "Input"],
Cell[70281, 1449, 13918, 243, 3190, "Input"],
Cell[84202, 1694, 361, 6, 170, "Input"],
Cell[84566, 1702, 694, 11, 270, "Input"],
Cell[85263, 1715, 432, 7, 190, "Input"],
Cell[85698, 1724, 891, 16, 250, "Input"],
Cell[86592, 1742, 293, 4, 150, "Input"],
Cell[86888, 1748, 232, 4, 150, "Input"],
Cell[87123, 1754, 700, 12, 190, "Input"],
Cell[87826, 1768, 660, 12, 190, "Input"],
Cell[88489, 1782, 734, 12, 210, "Input"],
Cell[89226, 1796, 269, 4, 150, "Input"]
}
]
*)



(*******************************************************************
End of Mathematica Notebook file.
*******************************************************************)

