0

This is in astronomy, but I think my question is probably very elementary - I'm not very experienced, I apologise.

I am plotting the relationship between the colour of a star-forming galaxy (y axis) with the redshift (x axis). The plot is a line that rises up from around 0 up to maybe 9, then decays again to about -2. The peak (~9 colour) is around 4 in terms of redshift, and I want to find the peak is more exactly. The redshift is given by quite a confusing function, and I can't figure out how to differentiate it or else I would just do that.

  • Could I maybe differentiate the complicated redshift (z) function? If so, how?

  • If not, how could I estimate a peak graphically/numerically?

Sorry for the very basic question and thank you very much in advance. My code is below.

import numpy as np
import matplotlib.pyplot as plt
import IGM
import scipy.integrate as integrate

SF = np.load('StarForming.npy')

lam = SF[0]

SED = SF[1]

filters = ['f435w','f606w','f814w','f105w','f125w','f140w','f160w']
filters_wl = {'f435w':0.435,'f606w':0.606,'f814w':0.814,'f105w':1.05,'f125w':1.25,'f140w':1.40,'f160w':1.60} # filter dictionary to give wavelengths of filters in microns

fT = {} # this is a dictionary

for f in filters:

data = np.loadtxt(f+'.txt').T

fT[f]= data    

fluxes = {}

for f in filters: fluxes[f] = [] # make empty list for each 

redshifts = np.arange(0.0,10.0,0.1) # redshifts going from 0 to 10

for z in redshifts:

    lamz = lam * (1. + z)
    obsSED = SED * IGM.madau(lamz, z)

    for f in filters:

        newT = np.interp(lamz,fT[f][0],fT[f][1]) # for each filter, refer back

        bb_flux = integrate.trapz((1./lamz)*obsSED*newT,x=lamz)/integrate.trapz((1./lamz)*newT,x=lamz) 
        # 1st bit integrates, 2nd bit divides by area under filter to normalise filter

        # loops over all z, for all z it creates a new SED, redshift wl grid    

        fluxes[f].append(bb_flux)


for f in filters: fluxes[f] = np.array(fluxes[f])

colour = -2.5*np.log10(fluxes['f435w']/fluxes['f606w'])

plt.plot(redshifts,colour)
plt.xlabel('Redshift')
plt.ylabel('Colour')
plt.show
albc
  • 23
  • 2
  • 6
  • For anyone who might want to try out your code, where does the IGM module come from? Also, would you please indent your code so that it works correctly. – Bill Bell Nov 10 '16 at 20:13
  • where is some sample data and its plots for visual verification (if anyone would want to try to code it)? – Spektre Nov 11 '16 at 10:08

1 Answers1

0

I do not have high enough reputation to comment, but this may solve your problem, so I guess its answer. Store all your y-coordinates in a list, then use the max(list) function to find the max. If you want an ordered pair, store your coordinates as (y,x) tuples and use max(list)

lst = [(3,2), (4,1), (1, 200)]
max(lst)

yields (4,1)

appills
  • 172
  • 1
  • 14
  • that is not precise enough... some interpolation reconstruction is needed to add to this (like intersecting the neighboring slopes similar to subpixel precision) the precision is curvature dependent so without data it is hard to say which interpolation to use (I mean order) I had in mind something like this: [Algorithm: how calculate INVERSE of bilinear interpolation?](http://stackoverflow.com/a/23103173/2521214). There are also another approaches like [Approximation search](http://stackoverflow.com/a/36163847/2521214) etc ... – Spektre Nov 11 '16 at 10:11