filter.c File Reference

Create radiometric filter functions. More...

#include "jurassic.h"

Go to the source code of this file.

Functions
double	ails (int apo, double opl, double dnu)
	Compute apodized instrument line shape. More...
int	main (int argc, char *argv[])

Detailed Description

Create radiometric filter functions.

Definition in file filter.c.

Function Documentation

ails()

double ails	(	int	apo,
		double	opl,
		double	dnu
	)

Compute apodized instrument line shape.

Definition at line 121 of file filter.c.

              {
 
  /* Auxiliary quantities... */
  const double a = 2 * M_PI * dnu * opl;
  const double a2 = a * a;
  const double a4 = a2 * a2;
  const double a6 = a4 * a2;
  const double a8 = a4 * a4;
 
  /* Sinc function... */
  if (apo == 0) {
    if (fabs(a) < 0.7)
      return 1 - a2 / 6 + a4 / 120 - a6 / 5040 + a8 / 362880;
    else
      return sin(a) / a;
  }
 
  /* Norton-Beer strong apodization... */
  else if (apo == 1) {
    double q0, q2, q4;
    if (fabs(a) < 0.7) {
      q0 = 1 - a2 / 6 + a4 / 120 - a6 / 5040 + a8 / 362880;
      q2 = 1 - a2 / 14 + a4 / 504 - a6 / 33264 + a8 / 3459456;
      q4 = 1 - a2 / 22 + a4 / 1144 - a6 / 102960 + a8 / 14002560;
    } else {
      const double sinca = sin(a) / a;
      const double cosa = cos(a);
      q0 = sinca;
      q2 = -15 / a2 * ((1 - 3 / a2) * sinca + (3 / a2) * cosa);
      q4 =
        945 / a4 * ((1 - 45 / a2 + 105 / a4) * sinca +
                    5 / a2 * (2 - 21 / a2) * cosa);
    }
    return 0.045335 * q0 + 0.554883 * q2 * 8. / 15. +
      0.399782 * q4 * 384. / 945.;
  }
 
  /* Error message.... */
  else
    ERRMSG("Unknown apodization!");
}

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 41 of file filter.c.

                {
 
  static ctl_t ctl;
 
  static double ff[NSHAPE], fnu[NSHAPE];
 
  double fsum = 0.0;
 
  int fn = 0;
 
  /* Write info... */
  if (argc < 3)
    ERRMSG("Give parameters: <ctl> <filter>");
 
  /* Read control parameters... */
  read_ctl(argc, argv, &ctl);
  const int type = (int) scan_ctl(argc, argv, "FILTER_TYPE", -1, "1", NULL);
  const double opd = scan_ctl(argc, argv, "FILTER_OPD", -1, "10.0", NULL);
  const double fwhm = scan_ctl(argc, argv, "FILTER_FWHM", -1, "1.0", NULL);
  const double center =
    scan_ctl(argc, argv, "FILTER_CENTER", -1, "1000.0", NULL);
  const double width = scan_ctl(argc, argv, "FILTER_WIDTH", -1, "2.1", NULL);
  const double samp = scan_ctl(argc, argv, "FILTER_SAMP", -1, "0.0005", NULL);
 
  /* Compute filter function... */
  for (double nu = center - 0.5 * width;
       nu <= center + 0.5 * width; nu += samp) {
 
    /* Set frequency... */
    fnu[fn] = nu;
 
    /* Boxcar... */
    if (type == 0)
      ff[fn] = (fabs(nu - center) <= 0.5 * fwhm ? 1.0 : 0.0);
 
    /* Triangle... */
    else if (type == 1) {
      ff[fn] = 1.0 - fabs(nu - center) / fwhm;
      ff[fn] = MAX(ff[fn], 0.0);
    }
 
    /* Gaussian... */
    else if (type == 2) {
      double sigma = fwhm / 2.355;
      ff[fn] = exp(-0.5 * POW2((nu - center) / sigma));
    }
 
    /* Sinc function... */
    else if (type == 3)
      ff[fn] = ails(0, opd, nu - center);
 
    /* Norton-Beer strong apodization... */
    else if (type == 4)
      ff[fn] = ails(1, opd, nu - center);
 
    /* Error message... */
    else
      ERRMSG("Filter function type unknown!");
 
    /* Count spectral grid points... */
    if ((++fn) >= NSHAPE)
      ERRMSG("Too many filter function data points!");
  }
 
  /* Normalize filter function... */
  for (int i = 0; i < fn; i++)
    fsum += ff[i];
  for (int i = 0; i < fn; i++)
    ff[i] /= (fsum * samp);
 
  /* Write to file... */
  write_shape(argv[2], fnu, ff, fn);
 
  return (EXIT_SUCCESS);
}

Here is the call graph for this function: