GASP_CoDeDi_C_code - david-macmahon/wiki_convert_test GitHub Wiki

/* ================================================================== */
/* Function to compute the dedispersion convolution phase function
   (Chirp function).
      DM         -> dispersion measure (pc/cc)
      BW         -> bandwidth (MHz)
      FCen       -> Centre frequency of the band (MHz)
      BRevSwitch -> Band Reversal switch (1 -> forward; -1 -> reversed)
      Chirp[]    -> Output chirp function array
      NPts       -> Number of points in Chirp[] (= NPts in FFT)

  While computing the phase value, it has been assumed that:
     charge of electron = 4.80325E-10 e.s.u.
     mass of electron   = 9.10956E-28 gm.
     velocity of light  = 2.997925E+10 cm/s
     parsec             = 3.0856E+18 cm.

            R. Ramachandran, Paul Demorest, 19-APR-2004, Berkeley     */
/* ================================================================== */

#include <stdio.h>
#include <math.h>
#include <fftw3.h>


int CalcChirp(double         DM,
              double         BW,
              double         Fcen,
              int            BRevSwitch,
              fftwf_complex  Chirp[],
              double         EarthZ4,
              int            NPts) {

  int     i, retval=0;
  float   *chr, *chi;
  double  Flow, dmconst, freq, frqfac, chw, pi, twopi, arg;
  float   taper;

  retval   = 0;
  pi       = 4.0*atan2(1.0,1.0);
  twopi    = (2.0 * pi);
  dmconst  = DM * twopi / ( 2.41e-10 * (1.0+EarthZ4/10000.0));  /* Earth correction */

  chw      = (BW / (float)NPts);
  Flow     = (Fcen - BW / 2.0) + (chw / 2.0);

  chr      = (float *)(Chirp);
  chi      = (float *)(Chirp) + 1;

  for (i=0; i<NPts; i++) {
    freq   = chw * ((double)i);
    if (i > NPts/2) { freq -= BW; }
    frqfac = (freq / Fcen);
    frqfac = (double)BRevSwitch*frqfac*frqfac
        /(Fcen+(double)BRevSwitch*freq);
    arg    = (dmconst * frqfac);
    taper = 1.0/sqrt(1.0 + pow((freq/(0.47*BW)),80));

    *chr   = (float)(cos(arg)*taper/(double)NPts); chr += 2;
    *chi   = -1.0*(float)(sin(arg)*taper/(double)NPts); chi += 2;
  }

  retval = 1;
  return retval;
}