chngres.c

Go to the documentation of this file.

/* 
Copyright 2021 Arizona Board of Regents on behalf of the University of Arizona.

For commercial uses, to obtain a license to sell and or sublicense copies of
the Software please contact the University of Arizona at Tech Launch Arizona:
lewish@tla.arizona.edu.

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction for non-commercial use educational and
research purposes, including without limitation the rights to use, copy,
modify, merge, publish, distribute, and to permit persons to whom the Software
is furnished to do so, subject to the following conditions:

The entire above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

1. Additional Required Provisions

1.1.  Arbitration. The parties agree that if a dispute arises between them
concerning this Agreement, the parties may be required to submit the matter to
arbitration pursuant to Arizona law.

1.2.  Applicable Law and Venue. This Agreement shall be interpreted pursuant
to the laws of the State of Arizona. Any arbitration or litigation between the
Parties shall be conducted in Pima County, ARIZONA, and LICENSEE hereby
submits to venue and jurisdiction in Pima County, ARIZONA.

1.3.  Non-Discrimination. The Parties agree to be bound by state and federal
laws and regulations governing equal opportunity and non-discrimination and
immigration.

1.4.   Appropriation of Funds. The Parties recognize that performance by
ARIZONA may depend upon appropriation of funds by the State Legislature of
ARIZONA. If the Legislature fails to appropriate the necessary funds, or if
ARIZONA’S appropriation is reduced during the fiscal year, ARIZONA may cancel
this Agreement without further duty or obligation. ARIZONA will notify
LICENSEE as soon as reasonably possible after it knows of the loss of funds.

1.5.   Conflict of Interest. This Agreement is subject to the provisions of
A.R.S. 38-511 and other conflict of interest regulations. Within three years
of the EFFECTIVE DATE, ARIZONA may cancel this Agreement if any person
significantly involved in initiating, negotiating, drafting, securing, or
creating this Agreement for or on behalf of ARIZONA becomes an employee or
consultant in any capacity of LICENSEE with respect to the subject matter of
this Agreement.

 */


#include "extern.h"

#include <gsl/gsl_math.h>

/*
 *-------------------------------------------------------------------------------
 *  @(#)chgres.f        5.2 09/10/98
 *-------------------------------------------------------------------------------
 *
 *  THIS SUBROUTINE SMOOTHS SPECTRAL LINE DATA FROM ITS INITIAL
 *  RESOLUTION TO A LOWER RESOLUTION SPECIFIED BY THE ADVERB
 *  'NEWRES'.  THE ROUTINE WORKS BY STEPPING A GAUSSIAN OF
 *  APPROPRIATE FWHM THROUGH THE EXISTING SPECTRUM AND PERFORMING
 *  A CONVOLUTION AT EACH POSITION.  THE ROUTINE ASSUMES THAT
 *  THERE ARE AN EVEN NUMBER OF CHANNELS IN THE INITIAL SPECTRUM
 *  AND THAT THE CENTER FREQENCY IN THE SCAN HEADER CORRESPONDS
 *  TO CHANNEL N/2 + 0.5, E.G., 64.5, 128.5, ETC.  EACH INDIVIDUAL
 *  CONVOLUTION IS PERFORMED OUT TO A SPECIFIED NUMBER OF STANDARD
 *  DEVIATIONS OF THE CONVOLVING GAUSSIAN, SAY 5 OR 6 SIGMA.  THE
 *  SMOOTHING IS TERMINATED WHEN THE CONVOLVING GAUSSIAN IS WITHIN
 *  A SPECIFIED NUMBER OF SIGMAS OF THE ENDS OF THE INITIAL SPECTRUM
 *  SO THAT THE END POINTS OF THE NEW SPECTRUM REPRESENT VALID 
 *  INFORMATION (2 OR 3 SIGMA IS A REASONABLE LIMIT).  THE REMAINING
 *  POINTS IN THE NEW SPECTRUM ARE SET EQUAL TO ZERO;  THE USER CAN
 *  DROP THESE POINTS FROM THE DISPLAYED SPECTRUM.
 *
 *  The center freq./reference channel need not be N/2 + 0.5, the
 *   reference channel is corrected in au2 after chgres returns
 *
 * ** WRITTEN BY P. R. JEWELL AND E. B. STOBIE  7 NOV 1985 **
 * ** LAST MODIFIED BY P.R.J. AND E.B.S. 9 NOV 1985 **
 *   Modified 8903 [RJM] See CHANGES.DOC
 *
 * ** Converted to 'C', April 12, 2016; t.w.folkers **
 *
 *  KEY TO VARIABLE NAMES:
 *      OLDRES = ORIGINAL FREQ. RESOLUTION OF THE SPECTRUM
 *      NCHAN = TOTAL NUMBER OF CHANNELS IN THE SPECTRUM
 *      NEWRES = DESIRED NEW FREQ. RESOLUTION OF THE SPECTRUM
 *             = FWHM of convolving gaussian
 *      CONRES = RESOLUTION OF THE CONVOLVING GAUSSIAN
 *      CONLIM = NUMBER OF SIGMAS OF THE CONVOLVING GAUSSIAN TO
 *               CARRY AN INDIVIDUAL CONVOLUTION OUT TO.
 *      CONEND = NUMBER OF SIGMAS OF THE CONVOLVING GAUSSIAN TO
 *               ALLOW AS A BORDER AT THE 2 ENDPOINTS OF THE SPEC.
 *
 */
int chgres(told, nchan, oldres, newres)
float *told, oldres, newres;
int nchan;
{
      int ncon, i, j, ichan, istart, istop;
      float tnew[8192];
      float conres, a2, rta2, conlim, conend,
            freqend, bord, conbord, cchan, chstep, chpos, xarg, xnorm,
            freqlc, chlim, xdist, arg, rinfinity = GSL_POSINF;

//      newres must be > oldres
      if(newres < oldres || min(newres, oldres) < 0.0) 
      {
        st_fail("CHNGRES: newres %f must be > oldres %f", newres, oldres);
        return(1);
      }

//      NCHAN must be even
      if((nchan%2) != 0)
      {
        st_fail("CHNGRES: nchan %d must be even", nchan);
        return(1);
      }

      if(nchan < 1 || nchan > 8196)
      {
        st_fail("CHNGRES: nchan %d out of range", nchan);
        return(1);
      }

      cchan = (float)nchan / 2.0 + 0.5;

/*
 *
 *  THE FWHM'S OF CONVOLVING GAUSSIANS ADD QUADRATICALLY:
 *      CONRES = SQRT( NEWRES*NEWRES - OLDRES*OLDRES )
 *      I don't understand why the below is preferred over the above - rwg.
 */
      conres = newres;
      a2     = 4.0 * log(2.0);
      rta2   = sqrt(a2);
      xarg   = rta2 * oldres / conres;
      conlim = 4.0;
      conend = 1.0;

/*
 *  A STANDARD DEVIATION SIGMA = 0.4248 * FWHM
 *                      0.4248 = 1 / sqrt(8 ln(2))
 *  CALCULATE THE NUMBER OF CONVOLUTIONS PERFORMED ON EITHER SIDE OF THE
 *    CENTER POINT
 */
      freqend = (nchan - cchan) * oldres;
      conbord = conend * 0.4248 * conres;
      ncon    = freqend / newres;
      freqlc  = ncon * newres;
      bord    = freqend - freqlc;

      if(bord < conbord)
      {
        ncon = ncon - 1;
      }
/*
 *
 *  IN THE LOOPS BELOW, THE OUTER LOOP STEPS THE CONVOLVING GAUSSIAN
 *  THROUGH THE SPECTRUM WITH STEPSIZES = NEWRES;  THE INNER LOOP
 *  PERFORMS THE INDIVIDUAL CONVOLUTIONS.  TWO SETS OF LOOPS ARE USED,
 *  THE FIRST WORKING FROM THE CENTER TO HIGHER CHANNELS AND THE
 *  SECOND FROM THE CENTER TO LOWER CHANNELS.
 *
 *  COMPUTE THE STEPSIZE, THE LIMIT OF THE CONVOLUTION, AND THE POSITION
 *  OF THE CENTER OF THE CONV. GAUSSIAN IN FRACTIONAL CHANNEL NUMBERS.
 *
 */
      chstep = newres / oldres;
      chlim  = conlim * 0.4248 * chstep;
      chpos  = cchan + 0.5 * chstep;

//  Initialize the new spectrum to zero
    bzero((char *)&tnew, sizeof(tnew));

/*
 *
 *  DO THE CONVOLUTION FOR THE UPPER HALF OF THE SPECTRUM:
 *
 */
      ichan = cchan;

      for(i=0;i<ncon;i++)
      {
        ichan = ichan + 1;
        istart = chpos - chlim + 0.5;
        istop = chpos + chlim + 0.5;

        if(istart < 1)
        {
          istart = 1;
        }

        if(istop > nchan)
        {
          istop = nchan;
        }

        xnorm = 0.0;

        for(j=istart;j<istop;j++)
        {
          xdist = chpos - j;
          arg = -xarg * xdist * xarg * xdist;

          if(gsl_finite(told[j]))
          {
            tnew[ichan] = told[j] * exp(arg) + tnew[ichan];
            xnorm += exp(arg);
          }
        }

        if(xnorm > 0.0)
        {
          tnew[ichan] /= xnorm;
        }
        else
        {
          tnew[ichan] = rinfinity;
        }

        chpos += chstep;
      }

/*
 *
 *  DO THE CONVOLUTION FOR THE LOWER HALF OF THE SPECTRUM:
 *
 */
      chpos = cchan - 0.5 * chstep;
      ichan = cchan + 1.0;

      for(i=0;i<ncon;i++)
      {
        ichan  = ichan - 1;
        istart = chpos - chlim + 0.5;
        istop  = chpos + chlim + 0.5;

        if(istart < 1)
        {
          istart = 1;
        }

        if(istop > nchan)
        {
          istop = nchan;
        }

        xnorm = 0.0;

        for(j=istart;j<istop;j++)
        {
          xdist = chpos - j;
          arg = -xarg * xdist * xarg * xdist;

          if(gsl_finite(told[j]))
          {
            tnew[ichan] = told[j] * exp(arg) + tnew[ichan];
            xnorm = xnorm + exp(arg);
          }
        }

        if(xnorm > 0.0)
        {
          tnew[ichan] /= xnorm;
        }
        else
        {
          tnew[ichan] = rinfinity;
        }

        chpos = chpos - chstep;
      }
/*
 *  
 *  REPLACE THE OLD SPECTRUM AND HEADER VALUES WITH THE NEW VALUES:
 *
 */
      for(i=0;i<nchan;i++)
      {
        if(i < cchan-ncon || i > cchan+ncon)
        {
          told[i] = rinfinity;
        }
        else
        {
          told[i] = tnew[i];
        }
      }

  return(0);
}