def fit_gauss(ph):
    for i in range(4):
		# bepaal de counts en de bingrootte in een range van 0 t/m 1000 ADC
        y, bins =  histogram(ph[i], bins=500, range=[0, 1000]) # loop for elke plaat
        x = (bins[:-1] + bins[1:])/2.  #let op neem bins om het midden te vinden van de bin
        
        if y[np.random.random_integers(200, 300)] != 0 : # ter controle of de data van de detector betrouwbaar is, aangezien plaat 4 503 geen goede data bevatte 3-01-2013
            # vanwege bartels (2012) kiezen we waar de ADC waarde tussen 150 en 410 zit voor de gauss fit, schat geeft de elementen waar tussen gefit moet worden.
            schat = np.where((x > 150) & (x < 600)) # geeft een tweeledig array terug    
            x1 = x[schat[0]] # geeft de waarde van de elementen gevonden met de determinatie hierboven in het eerste element
            #bovenstaande kan ook met x1 = fit_xa[np.where(x>150) & (x< 420)]
            y1 = y[schat] # bepaling van de count in de meting 
            max_min = ndimage.extrema(y1)
            max_y = ndimage.extrema(y1)[1] #maximale count y waarde piek van de gauss!)
        
            min_x = max_min[0]  #  de laagste waarde van de count
            max_x = max_min[1]  # hoogste waarde van de count -> a waarde in de gauss fit
            b_temp = max_min[3]
            b_place = b_temp[0] # b waarde
            b = x1[b_place] # de b-waarde van de gauss curve
        
            bound1 = max_x - (max_x - min_x)*0.75 
        
            if (max_x- min_x) <= 50:
                bound2 = max_x + (max_x - min_x)*1.5
            else:
                bound2 = max_x + (max_x - min_x)
        
            x2 = x1.compress((bound1 <= x1) & (x1 < bound2))
            y2 = y1.compress((bound1 <= x1) & (x1 < bound2))
        
              
            #popt, pcov = curve_fit(func1, x1, y1, [200, 250, 50] ) # werkende fit met een handmatige guess kleine verschillen. orde van 0.2 
            popt, pcov = curve_fit(func1, x1, y1, [max_x, b, 20] ) # fit met guess gebruikt werkt ook
        
            pylab.plot(func1(range(0,600), *popt), '--')
            peak = popt[1]
        
            if i == 0:
                MIP1 = peak
                MPV1.append(MIP1)  #bij herhalen  van de loop oor ander station append deze regel op de juiste manier?
            elif i == 1:
                MIP2 = peak
                MPV2.append(MIP2)
            elif i == 2:
                MIP3 = peak
                MPV3.append(MIP3)
            elif i == 3:
                MIP4 = peak
                MPV4.append(MIP4)
        
            print 'The MPV of detector',i + 1, 'lies at', peak, 'ADC'
        
        else:
            
            print 'The data of the detector ',i + 1, 'could not be fitted to a gauss curve'

    def new_image(self, data):
        """
        Update the window with the new image (the window is resize to have the image
        at ratio 1:1)
        data (numpy.ndarray): an 2D array containing the image (can be 3D if in RGB)
        """
        if data.ndim == 3 and 3 in data.shape: # RGB
            rgb = img.ensureYXC(data)
        elif numpy.prod(data.shape) == data.shape[-1]: # 1D image => bar plot
            # TODO: add "(plot)" to the window title
            # Create a simple bar plot of X x 400 px
            lenx = data.shape[-1]
            leny = 400
            maxy = data.max()
            logging.info("Plot data max = %s", maxy)
            rgb = numpy.zeros((leny, lenx, 3), dtype=numpy.uint8)
            for i, v in numpy.ndenumerate(data):
                h = leny - int((v * leny) / maxy)
                rgb[h:-1, i[-1], :] = 255
        else: # Greyscale (hopefully)
            mn, mx, mnp, mxp = ndimage.extrema(data)
            logging.info("Image data from %s to %s", mn, mx)
            rgb = img.DataArray2RGB(data) # auto brightness/contrast

        self.app.img = NDImage2wxImage(rgb)
        wx.CallAfter(self.app.update_view)

def normalized_per_object(image, labels):
    """Normalize the intensities of each object to the [0, 1] range."""
    nobjects = labels.max()
    objects = np.arange(nobjects + 1)
    lmin, lmax = scind.extrema(image, labels, objects)[:2]
    # Divisor is the object's max - min, or 1 if they are the same.
    divisor = np.ones((nobjects + 1,))
    divisor[lmax > lmin] = (lmax - lmin)[lmax > lmin]
    return (image - lmin[labels]) / divisor[labels]

 def new_image(self, data):
     """
     Update the window with the new image (the window is resize to have the image
     at ratio 1:1)
     data (numpy.ndarray): an 2D array containing the image (can be 3D if in RGB)
     """
     mn, mx, mnp, mxp = ndimage.extrema(data)
     logging.info("Image data from %s to %s", mn, mx)
     rgb = DataArray2RGB(data) # auto brightness/contrast
     self.app.img = NDImage2wxImage(rgb)
     wx.CallAfter(self.app.update_view)

 def f(array, il):
     i = il[0]
     # This function returns (mini, maxi), the indexes of the max
     # and min of the array.  They return the *lowest* possible
     # such indexes, thus the max(0, ...) below.
     min_, max_, mini, maxi = extrema(array)
     #print array, (min_, max_, mini, maxi), mini[0]+i-halfwidth, maxi[0]+i-halfwidth
     minima[max(0, mini[0]+i-halfwidth)] += 1
     maxima[       maxi[0]+i-halfwidth ] += 1
     il[0] += 1
     return 0

def test_extrema01():
    labels = np.array([1, 0], bool)
    for type in types:
        input = np.array([[1, 2], [3, 4]], type)
        output1 = ndimage.extrema(input, labels=labels)
        output2 = ndimage.minimum(input, labels=labels)
        output3 = ndimage.maximum(input, labels=labels)
        output4 = ndimage.minimum_position(input,
                                                     labels=labels)
        output5 = ndimage.maximum_position(input,
                                                     labels=labels)
        assert_equal(output1, (output2, output3, output4, output5))

def scipy_features(image):
    histo = ndimage.measurements.histogram(image, 0, 255, 10)
    mini, maxi, minpos, maxpos = ndimage.extrema(image)
    return numpy.concatenate((
        histo,
        minpos,
        maxpos,
        ndimage.measurements.center_of_mass(image),
        [
            mini,
            maxi,
        ]
    ))

def grad_graph(img, interactive=True):
    """
    gradient graph: use gradient to define different objects and then segment the image to those objects
    :param img: image to process
    :return: fragmentation
    """

    blur_radius = 2
    imgf = ndimage.gaussian_filter(img, blur_radius)


    scharr2 = np.array([[ -3, 0,  3],
                       [-10, 0, 10],
                       [ -3, 0,  3]])
    grad = [convolve2d(imgf[:, :, i], scharr2, boundary='symm', mode='same') for i in np.arange(3)]
    grad = np.sum(np.abs(grad), 0, dtype=np.uint8)
    imgf[ndimage.gaussian_filter(grad/np.mean(grad), 3) < 1] = 0
    threshold = np.mean(imgf)
    mask = np.any(imgf > threshold, -1)
    labeled, nr_objects = ndimage.label(mask)
    # labeled: a width X height matrix, in which each pixel of object i is given index i

    sizes = ndimage.sum(mask, labeled, range(nr_objects + 1))
    #plt.title('objects are colored')
    #plt.imshow(labeled)
    #plt.show()
    labeled[(sizes < 1000)[labeled]] = 0  # filter noise
    label_im = np.searchsorted(np.unique(labeled), labeled)
    # TODO: can we do better with widths and heights (Eran: I think it is the way to calc width/height)
    widths_heights = np.abs(np.diff((np.array(ndimage.extrema(mask, label_im, range(np.max(label_im) + 1))[2:])), axis=0)[0, :, :])
    label_im[(widths_heights[:, 0] < 5) | (widths_heights[:, 0] > 800) | (widths_heights[:, 1] > 300)] = 0
    label_im = np.searchsorted(np.unique(label_im), label_im)

    #plt.hist(widths_heights[:, 0],100) width histogram
    #plt.show()

    # now given the objects return their positions (frags)
    frags = ndimage.find_objects(label_im, max_label=np.max(label_im))
    if interactive:
        # show example of one of the words
        img2 = np.array(img, copy=True)
        img2[label_im == 0] = 0
        plt.imshow(img2)
        plt.show()

        #
        fragment_to_show = len(frags)//2
        plt.imshow(img[frags[fragment_to_show]])
        plt.show()

    return frags

    def new_image(self, data):
        """
        Update the window with the new image (the window is resize to have the image
        at ratio 1:1)
        data (numpy.ndarray): an 2D array containing the image (can be 3D if in RGB)
        """
        if data.ndim == 3 and 3 in data.shape: # RGB
            rgb = img.ensureYXC(data)
        else: # Greyscale (hopefully)
            mn, mx, mnp, mxp = ndimage.extrema(data)
            logging.info("Image data from %s to %s", mn, mx)
            rgb = img.DataArray2RGB(data) # auto brightness/contrast

        self.app.img = NDImage2wxImage(rgb)
        wx.CallAfter(self.app.update_view)

def test_extrema02():
    "extrema 2"
    labels = np.array([1, 2])
    for type in types:
        input = np.array([[1, 2], [3, 4]], type)
        output1 = ndimage.extrema(input, labels=labels,
                                            index=2)
        output2 = ndimage.minimum(input, labels=labels,
                                            index=2)
        output3 = ndimage.maximum(input, labels=labels,
                                            index=2)
        output4 = ndimage.minimum_position(input,
                                            labels=labels, index=2)
        output5 = ndimage.maximum_position(input,
                                            labels=labels, index=2)
        assert_equal(output1, (output2, output3, output4, output5))

def test_extrema04():
    labels = [1, 2, 0, 4]
    for type in types:
        input = np.array([[5, 4, 2, 5],
                                [3, 7, 8, 2],
                                [1, 5, 1, 1]], type)
        output1 = ndimage.extrema(input, labels, [1, 2])
        output2 = ndimage.minimum(input, labels, [1, 2])
        output3 = ndimage.maximum(input, labels, [1, 2])
        output4 = ndimage.minimum_position(input, labels,
                                                     [1, 2])
        output5 = ndimage.maximum_position(input, labels,
                                                     [1, 2])
        assert_array_almost_equal(output1[0], output2)
        assert_array_almost_equal(output1[1], output3)
        assert_array_almost_equal(output1[2], output4)
        assert_array_almost_equal(output1[3], output5)

    def mark_peaks_of_trace(self, absolute=None, relative=None):
        '''Absolute means that the peak-finding is based on the raw value.  
           Relative is a fraction of the maximum. '''
        # Calculate the cutoff
        # The logic is inflated for sanity check 
        if (absolute is not None) and (relative is None):
            ekg_cutoff = absolute
        elif (relative is not None) and (absolute is None):
            ekg_cutoff = self.global_ekg_max * relative
        else:
            # Probably not necessary, as ekg_cutoff would be undefined
            raise RuntimeError("Must sepcify ekg cutoff in one way!")
        
        # Stuff to store internally. 
        self.spike_times = []   # The times in EKG ms offset
        self.spike_number = []  # The number of spikes
        self.spike_heights = [] # The amplitudes

        # For each EKG trace
        for n, tr in enumerate(self.traces):
            # Find the area above the cutoff
            screen = (tr > ekg_cutoff)

            # Label each and record the number of labels
            labels, spike_count = ndimage.label(screen)
            
            # Find the extrema; locations and amplitudes
            indexs_of_spikes = range(1, spike_count + 1)

            # This conditional should not be necessary.  
            # TODO: scipy bug-tracker is down . . .  report this one
            if indexs_of_spikes != []:
                ext = ndimage.extrema(tr, labels=labels, index=indexs_of_spikes)
            else:
                ext = ([],[],[],[])
            
            # Break out the results of the extrema, ignore min-related values
            min_vals, max_vals, min_locs, max_locs = ext

            # Store the values on a per-experiment basis
            # The arrays appended _could_ be empty, but will be an object . . . 
            # so the offsets will still be correct.
            self.spike_times.append( array(max_locs).flatten() )
            self.spike_number.append( spike_count )
            self.spike_heights.append( array(max_vals).flatten() ) 

def test_extrema03():
    labels = np.array([[1, 2], [2, 3]])
    for type in types:
        input = np.array([[1, 2], [3, 4]], type)
        output1 = ndimage.extrema(input, labels=labels,
                                            index=[2, 3, 8])
        output2 = ndimage.minimum(input, labels=labels,
                                            index=[2, 3, 8])
        output3 = ndimage.maximum(input, labels=labels,
                                            index=[2, 3, 8])
        output4 = ndimage.minimum_position(input,
                                    labels=labels, index=[2, 3, 8])
        output5 = ndimage.maximum_position(input,
                                    labels=labels, index=[2, 3, 8])
        assert_array_almost_equal(output1[0], output2)
        assert_array_almost_equal(output1[1], output3)
        assert_array_almost_equal(output1[2], output4)
        assert_array_almost_equal(output1[3], output5)

 def new_image(self, data):
     """
     Update the window with the new image (the window is resize to have the image
     at ratio 1:1)
     data (numpy.ndarray): an 2D array containing the image (can be 3D if in RGB)
     """
     if data.ndim == 3 and 3 in data.shape:  # RGB
         rgb = img.ensureYXC(data)
     elif numpy.prod(data.shape) == data.shape[-1]:  # 1D image => bar plot
         # TODO: add "(plot)" to the window title
         # Create a simple bar plot of X x 400 px
         lenx = data.shape[-1]
         if lenx > MAX_WIDTH:
             binning = lenx // MAX_WIDTH
             data = data[..., 0::binning]
             logging.debug("Compressed data from %d to %d elements", lenx, data.shape[-1])
             lenx = data.shape[-1]
         leny = 400
         miny = min(0, data.min())
         maxy = data.max()
         diffy = maxy - miny
         if diffy == 0:
             diffy = 1
         logging.info("Plot data from %s to %s", miny, maxy)
         rgb = numpy.zeros((leny, lenx, 3), dtype=numpy.uint8)
         for i, v in numpy.ndenumerate(data):
             # TODO: have the base at 0, instead of miny, so that negative values are columns going down
             h = leny - int(((v - miny) * leny) / diffy)
             rgb[h:-1, i[-1], :] = 255
     else:  # Greyscale (hopefully)
         mn, mx, mnp, mxp = ndimage.extrema(data)
         logging.info("Image data from %s to %s", mn, mx)
         rgb = img.DataArray2RGB(data)  # auto brightness/contrast
     self.app.spots, self.app.translation, self.app.scaling, self.app.rotation = FindGridSpots(data, self.gridsize)
     self.app.img = NDImage2wxImage(rgb)
     wx.CallAfter(self.app.update_view)

#.........这里部分代码省略.........
                order = np.zeros(distance_matrix.shape, int)
            else:
                order = np.lexsort([distance_matrix])
            first_neighbor = 1 if self.neighbors_are_objects else 0
            first_object_index = order[:, first_neighbor]
            first_x_vector = ncenters[first_object_index,1] - ocenters[:,1]
            first_y_vector = ncenters[first_object_index,0] - ocenters[:,0]
            if nneighbors > first_neighbor+1:
                second_object_index = order[:, first_neighbor + 1]
                second_x_vector = ncenters[second_object_index,1] - ocenters[:,1]
                second_y_vector = ncenters[second_object_index,0] - ocenters[:,0]
                v1 = np.array((first_x_vector,first_y_vector))
                v2 = np.array((second_x_vector,second_y_vector))
                #
                # Project the unit vector v1 against the unit vector v2
                #
                dot = (np.sum(v1*v2,0) / 
                       np.sqrt(np.sum(v1**2,0)*np.sum(v2**2,0)))
                angle = np.arccos(dot) * 180. / np.pi
            
            # Make the structuring element for dilation
            strel = strel_disk(distance)
            #
            # A little bigger one to enter into the border with a structure
            # that mimics the one used to create the outline
            #
            strel_touching = strel_disk(distance + .5)
            #
            # Get the extents for each object and calculate the patch
            # that excises the part of the image that is "distance"
            # away
            i,j = np.mgrid[0:labels.shape[0],0:labels.shape[1]]
            min_i, max_i, min_i_pos, max_i_pos =\
                scind.extrema(i,labels,object_indexes)
            min_j, max_j, min_j_pos, max_j_pos =\
                scind.extrema(j,labels,object_indexes)
            min_i = np.maximum(fix(min_i)-distance,0).astype(int)
            max_i = np.minimum(fix(max_i)+distance+1,labels.shape[0]).astype(int)
            min_j = np.maximum(fix(min_j)-distance,0).astype(int)
            max_j = np.minimum(fix(max_j)+distance+1,labels.shape[1]).astype(int)
            #
            # Loop over all objects
            # Calculate which ones overlap "index"
            # Calculate how much overlap there is of others to "index"
            #
            for object_number in object_numbers:
                if object_number == 0:
                    #
                    # No corresponding object in small-removed. This means
                    # that the object has no pixels, e.g. not renumbered.
                    #
                    continue
                index = object_number - 1
                patch = labels[min_i[index]:max_i[index],
                               min_j[index]:max_j[index]]
                npatch = neighbor_labels[min_i[index]:max_i[index],
                                         min_j[index]:max_j[index]]
                #
                # Find the neighbors
                #
                patch_mask = patch==(index+1)
                extended = scind.binary_dilation(patch_mask,strel)
                neighbors = np.unique(npatch[extended])
                neighbors = neighbors[neighbors != 0]
                if self.neighbors_are_objects:
                    neighbors = neighbors[neighbors != object_number]

def find_MPV(fit, conv=1):
    """Determine the ADC value of the MPV for each detector

    For pulseheights use conv = 1
    For pulseintegral values use conv = 10

    Returns a tuple of MPV values

    """
    MPV1, MPV2, MPV3, MPV4 = [], [], [], []

    for i in range(4): ##4 DETECTOREN
        try:
            y, bins, patches = plt.hist(fit[i], bins=linspace(0, conv * 2000 / 0.57, STEPSIZE), label="Data")

            x = bins[:-1] + .5 * (bins[1] - bins[0])

            # A first approximation of the position of the MPV and the minimum to the left is made. Used for fit bounds and initial guess of the parameters.
            guess_max_x = ndimage.extrema(y.compress((conv * 150 <= x) & (x < conv * 410)))[3][0] + len(y.compress(x < conv * 150))
            guess_max_x = x[guess_max_x]
            guess_min_x = ndimage.extrema(y.compress((conv * 53 <= x) & (x < guess_max_x)))[2][0] + len(y.compress(x < conv * 53))
            guess_min_x = max([conv * 120, x[guess_min_x]])
            guess_max_y = ndimage.extrema(y.compress((conv * 150 <= x) & (x < conv * 410)))[1]
            guess_min_y = ndimage.extrema(y.compress((conv * 53 <= x) & (x < guess_max_x)))[0]

            g0 = guess_max_y
            g1 = guess_max_x

            # The fit range. Since the right side of the ph histogram is most similar to a Gauss function, we do not want to set bound 2 too small.
            bound1 = guess_min_x + (guess_max_x - guess_min_x) / 4

            if (guess_max_x - guess_min_x) <= conv * 50:
                bound2 = guess_max_x + (guess_max_x - guess_min_x) * 1.5
            else:
                bound2 = guess_max_x + (guess_max_x - guess_min_x)

            # Fit a Gaussian to the peak.
            f = lambda x, a, b, c: a * exp(-((x - b) ** 2) / c)
            x2 = x.compress((bound1 <= x) & (x < bound2))
            y2 = y.compress((bound1 <= x) & (x < bound2))
            popt, pcov = optimize.curve_fit(f, x2, y2, (g0, g1, (guess_max_x - guess_min_x) / 2), sigma=sqrt(y2))

            # Find the x value of the peak:
            peak = popt[1]
            print ('The MPV of the pulseheight of detector %d lies at %d ADC' %
                   (i + 1, peak))

        # Safety nets: prevent the code from terminating if an error occurs
        except IndexError:
            peak = nan
            print 'Apparently there is no data'
        except RuntimeError:
            peak = nan
            print 'Unable to make a good fit'
        except TypeError:
            peak = nan
            print 'Not enough data to make a proper fit' # When it uses less than three bins to make the fit you get this error
        except ValueError:
            peak = nan
            print 'Value Error Occured'

        if i == 0:
            MIP1 = peak
        elif i == 1:
            MIP2 = peak
        elif i == 2:
            MIP3 = peak
        elif i == 3:
            MIP4 = peak

    mpv = (MIP1, MIP2, MIP3, MIP4)
    return mpv

 def measure_worms(self, workspace, labels, nworms, width):
     m = workspace.measurements
     assert isinstance(m, cpmeas.Measurements)
     object_name = self.straightened_objects_name.value
     nbins_vertical = self.number_of_segments.value
     nbins_horizontal = self.number_of_stripes.value
     params = self.read_params(workspace)
     if nworms == 0:
         # # # # # # # # # # # # # # # # # # # # # #
         #
         # Record measurements if no worms
         #
         # # # # # # # # # # # # # # # # # # # # # #
         for ftr in (FTR_MEAN_INTENSITY, FTR_STD_INTENSITY):
             for group in self.images:
                 image_name = group.straightened_image_name.value
                 if nbins_vertical > 1:
                     for b in range(nbins_vertical):
                         measurement = "_".join(
                             (C_WORM, ftr, image_name, 
                              self.get_scale_name(None, b)))
                         m.add_measurement(object_name, measurement,
                                           np.zeros((0)))
                 if nbins_horizontal > 1:
                     for b in range(nbins_horizontal):
                         measurement = "_".join(
                             (C_WORM, ftr, image_name, 
                              self.get_scale_name(b, None)))
                         m.add_measurement(object_name, measurement,
                                           np.zeros((0)))
                     if nbins_vertical > 1:
                         for v in range(nbins_vertical):
                             for h in range(nbins_horizontal):
                                 measurement = "_".join(
                                     (C_WORM, ftr, image_name, 
                                      self.get_scale_name(h, v)))
                                 m.add_measurement(object_name, measurement,
                                                   np.zeros((0)))
                                 
     else:
         #
         # Find the minimum and maximum i coordinate of each worm
         #
         object_set = workspace.object_set
         assert isinstance(object_set, cpo.ObjectSet)
         orig_objects = object_set.get_objects(self.objects_name.value)
 
         i,j = np.mgrid[0:labels.shape[0], 0:labels.shape[1]]
         min_i, max_i, _, _ = extrema(i, labels, orig_objects.indices)
         min_i = np.hstack(([0], min_i))
         max_i = np.hstack(([labels.shape[0]], max_i)) + 1
         heights = max_i - min_i
         
         # # # # # # # # # # # # # # # # #
         #
         # Create up to 3 spaces which represent the gridding
         # of the worm and create a coordinate mapping into
         # this gridding for each straightened worm
         #
         # # # # # # # # # # # # # # # # #
         griddings = []
         if nbins_vertical > 1:
             scales = np.array([self.get_scale_name(None, b)
                                for b in range(nbins_vertical)])
             scales.shape = (nbins_vertical, 1)
             griddings += [(nbins_vertical, 1, scales)]
         if nbins_horizontal > 1:
             scales = np.array([self.get_scale_name(b, None)
                                for b in range(nbins_horizontal)])
             scales.shape = (1, nbins_horizontal)
             griddings += [(1, nbins_horizontal, scales)]
             if nbins_vertical > 1:
                 scales = np.array([
                     [self.get_scale_name(h,v) for h in range(nbins_horizontal)]
                     for v in range(nbins_vertical)])
                 griddings += [(nbins_vertical, nbins_horizontal, scales)]
         
         for i_dim, j_dim, scales in griddings:
             # # # # # # # # # # # # # # # # # # # # # #
             #
             # Start out mapping every point to a 1x1 space
             #
             # # # # # # # # # # # # # # # # # # # # # #
             labels1 = labels.copy()
             i,j = np.mgrid[0:labels.shape[0], 0:labels.shape[1]]
             i_frac = (i - min_i[labels]).astype(float) / heights[labels]
             i_frac_end = i_frac + 1.0 / heights[labels].astype(float)
             i_radius_frac = (i - min_i[labels]).astype(float) / (heights[labels] - 1)
             labels1[(i_frac >= 1) | (i_frac_end <= 0)] = 0
             # # # # # # # # # # # # # # # # # # # # # #
             #
             # Map the horizontal onto the grid.
             #
             # # # # # # # # # # # # # # # # # # # # # #
             radii = np.array(params.radii_from_training)
             #
             # For each pixel in the image, find the center of its worm
             # in the j direction (the width)
             #
             j_center = int(width / 2) + width * (labels-1)
#.........这里部分代码省略.........

"""
Find at what time the temperature suddenly took off
"""

import sys, glob
import pyfits
from scipy.ndimage import extrema

runid = sys.argv[1]
pattern = "%s-te????.fits" % (runid)
filelist = glob.glob(pattern)

for file in filelist:
    Tarray = pyfits.open(file)[0].data
    Tmin, Tmax, Tminloc, Tmaxloc = extrema(Tarray)
    # switch to fortran array indexing
    Tmaxloc = [ i+1 for i in Tmaxloc[::-1] ]
    print "%s: max T = %.2e at %s" % (file, Tmax, Tmaxloc)