next up previous
Next: Updating magmax, filter and Up: The C code Previous: The IDL function CALL_EXTERNAL

  
Description of the functions in magnitude.so

The functions in the file magnitude.so are:

- magnitude
- readtable
- integrate
- linterpolate
- hunt

magnitude: The function that is called from the IDL function magmax is magnitude. In this function, first of all the parameters are converted as described above. Then the flag spec is set to 1 and readtable is called with the path of the SN spectrum. The flag spec tells readtable if a spectrum (1) or a filter (0) is to be read in (see 2.5.2, readtable below). The return value of readtable is an error flag, the arrays xspec0 and yspec0 contain the wavelength (in Å) an the flux (in photons), with an array offset of one. maxspec0 contains the length of these arrays. After reading the spectrum in, the error flag is tested. Then the spectrum is calibrated to the value of the parameter mb as magnitude in B band. For further reference on this calibration, see 2.6 below. The next step is reading in the filter, so the flag spec is set to 0. The variable fil contains an integer value indicating which filter to use (it is set by the calling IDL function). The correspondence is the following:

value of fil corresponding filter
1 U
2 B
3 V
4 R
5 I
6 J
7 H
8 K (short)

For more information on the filters, see section 2.4 above. After the filter has been read in, the error flag is checked again. The contents of the array x is the wavelength of the filter, y is the value at this wavelength and max contains the length of these two arrays. The normalisation of the filter is applied later (see 2.5.2, integrate below and also 2.6). Note that the J, H and Ks band have to be converted from nm to Å. The next step is to reserve memory for a copy of the spectrum. The original spectrum will be kept in memory and the copy will be used for working on the spectrum. This is necessary because the spectrum has to be redshifted. But if several redshifts are to be processed, the spectrum would have to be read in again for processing the next redshift. In order to save runtime, the original spectrum is kept in memory and only a copy of it will be redshifted.
The following for-loop contains the actual calculation. The for-loop makes sure that this calculation is done for every z value passed to magnitude. At the beginning of the for-loop, the original spectrum is copied to a working copy (see above). Then the copied spectrum is redshifted. After this, the integration routine integrate is called. It returns an error flag which is then checked. The result of the integration is stored in the array result. It contains the magnitude of the redshifted SN spectrum in the specified filter. For more information on integrate, see 2.5.2, integrate below. After the error flag handed back by integrate is checked, the distance modulus is added to the magnitude. Then the for-loop is entered again. After all z values have been processed by the for-loop, the memory is deallocated and magnitude returns to the IDL function magmax.

readtable:  This function is used to read in the spectrum and the filter and get the memory needed for this. First, the file specified by the string path is attempted to open. If this failes, readtable returns. If not, the first line is read. The data is assumed to be arranged in two columns in the file, the first one containing the wavelength, the second one containing the flux. The columns have to be separated by white space characters and nothing else is allowed within the file. After the first line has been read, the result of fscanf is checked. If it is 2, this means that 2 arguments have successfully been read and converted to the type indicated by the format string, i.e. double precision floating point. If the returned value wasn't 2, an error is produced, the file is closed and readtable returns to magnitude. If it was 2, readtable allocates memory for two arrays (the x and y values of the read table). At first, memory for only two double float numbers is allocated. x[0] and y[0] aren't used, so the first actual stored value will be x[1] and y[1]. After the memory allocation (and checking for its success) the values from the first read line are stored. If the flag spec indicates that a spectrum is read, the flux is converted to photons before it is stored. If spec indicates a filter, the read values are stored without changes. After this, a while loop is entered for reading in the rest of the file. The exit condition of this while loop is flag $\neq$ 1. In the loop, the next line is read and the result is checked if end-of-file is reached. If the end of the file is reached, flag is set to 0 so that the while loop exits. If it is not end of file, the read and conversion results are checked in the same way as described above. If the conversion was ok, counter is increased and realloc is called to get memory for one more double float, appended to the array x (or y). The returned pointer is checked for success. If the allocation was successful, the old array base pointer x (or y) is set to the new array. After both arrays have been modified this way, the read values are stored in the same way as before. After the entire file has been read in by the while loop, the length of the arrays is written into one of magnitude's variables, the file is closed and readtable returns to magnitude.

integrate:  This function performs the actual calculation. First, memory is allocated for a copy of the filter. If the allocation was successful, the spectrum is copied to the new arrays and the filter undergoes the same procedure. This is done because the interpolation routine used for the integration will change the arrays. The next step is checking if the filter doesn't begin at a lower wavelength than the spectrum or end at a higher wavelength than the spectrum. If so, integrate frees the allocated memory and returns to magnitude with an error. If the filter matches the spectrum, linterpolate is called to linearly interpolate the filter. The arguments are: array with x values of the function, array with y values for the function, length of the arrays, starting point of the interpolation, end point of the interpolation. For more information on this function, see 2.5.2, linterpolate below. The function linterpolate modifies the two arrays passed to it as well as their length. After the call to linterpolate, the arrays x and y contain the interpolated data beginning with the starting point specified before and ending with the ending point specified before. The length of the arrays has been adjusted to this. After the call to linterpolate, integrate checks the error flag handed back by linterpolate. Then the spectrum is interpolated in the same way. Note that the beginning / end point of the spectrum interpolation are identical with the beginning / end point of the filter. This is done to save runtime. The integration of the filter with the spectrum would cut out anything in the spectrum that is not matched by the filter. With the approach described above, the spectrum is truncated to the part which matches the filter and the rest is dropped. So only those values which really take effect in the integration are kept. Also, another feature of linterpolate is to interpolate spectrum and filter with the same stepsize in wavelength. So, as beginning / end point of spectrum and filter are identical, the integration reduces to just summing up the values in the filter * the (due to the identical stepsize, now corresponding) values inthe spectrum * the stepsize. This is the next step in function integrate. After the integration, the normalisation of the filters is applied. The data is stored in a struct and has been calculated with the program normfilter (see section 2.6 for more information on this). After all this, integrate returns to magnitude.

linterpolate:  This function is used by integrate to linearly interpolate the data for the integration. The parameters start and end are used as beginning / end point of the interpolation, l and f are the arrays containing wavelength and flux and max is their length. First of all, memory for the interpolated data is allocated. Again, just memory for two double floats is allocated, the first one isn't used, so the values in the arrays will start with *l[1], *f[1]. Of course, the success of the allocation is checked and linterpolate returns if this wasn't successful.
After this, a for-loop containing the actual interpolation is entered. This for-loop produces the desired interpolation points in wavelength, called x, running from start to end with a stepsize of dx. The stepsize dx is defined as a const double at the beginning of the source code because it is also used in integrate. In addition to this, a variable newmax is initialised. It will be used to store the length of the new, interpolated arrays. The increment of newmax, however, is placed into the loop itself. The reason for this will be explained later. In the loop, first the x value is checked if it is the beginning / end point of the wavelength array. If yes, the corresponding y value (to x) is just copied from the arrays instead of interpolating it. If x is not the beginning/ end of the wavelength, this means it is somewhere within the array, so interpolation makes sense. The first step in the interpolation is to search the position of x (or the position of the next lowest entry) in the wavelength array. This work is done by the function hunt. It was taken from the book ``Numerical recipes in C'', so for documentation on the function, check the book [1]. The function hunt needs a guess value for the position of x in the array, the mid of the array is not a bad start. The position is returned in pos. After that, the y value for x is linearly interpolated, y = m x + b. As the zeropoint of the wavelength is kind of shifted to *l[pos], x in that formula has to be the offset from *l[pos] to x and b is the flux value at pos. After the corresponding y value to x has been found, x and y have to be stored. If newmax is 1, there's still enough free memory in the allocated arrays newl and newf, so x and y can be stored there directly. Otherwise, new memory is allocated, it is appended to the existing arrays. If the increment of newmax would have been placed in the arguments of the for-statement, the memory management couldn't work like this. The success of the allocation is checked and if it was successful, the old array base pointers are used again to point to the new, enlarged arrays. The new memory is then used to store the interpolated values. After the whole interpolation is done, the memory containing the old wavelength and flux data is deallocated and the old pointers are used to store the base pointers of the interpolated data. Then linterpolate returns to integrate.

hunt: The function hunt is taken from the book ``Numerical recipes in C''. For more information on the function, check the book [1].


next up previous
Next: Updating magmax, filter and Up: The C code Previous: The IDL function CALL_EXTERNAL
Peter E. Nugent Jr.
1998-10-01