view release on metacpan or  search on metacpan

ast/src/fitschan.c  view on Meta::CPAN

   double **ptr1;
   double **ptr2;
   double *p;
   double *s;
   double *xl;
   double c;
   double d;
   double delta;
   double in_lbnd;
   double in_ubnd;
   double lbnd_out;
   double m;
   double p0;
   double pv;
   double sn;
   double sp;
   double sps;
   double ss2;
   double ss;
   double sv;
   double tol;
   double ubnd_out;
   int *ins;
   int boxok;
   int i;
   int j;
   int nin;
   int nout;
   int oldrep;
   int ret;

/* Initialise */
   ret = 0;
   if( coord_in ) *coord_in = -1;

/* Check inherited status */
   if( !astOK ) return ret;

/* Attempt to split off the required output (in case any of the other
   outputs are associated with Mappings that do not have an inverse). */
   astInvert( smap );
   ins = astMapSplit( smap, 1, &coord_out, &map );
   astInvert( smap );

/* If successful, check that the supplied output is fed by only one input.
   At the moment "map" goes form wcs to pixel, so we check its Nout
   attribute. If required, return the index of the input that feeds the
   requested output. */
   if( ins ) {
      if( astGetNout( map ) == 1 ) {
         if( coord_in ) *coord_in = ins[ 0 ];

/* If so, invert the map so that it goes from pixel to wcs. Here on we
   use "map" in place of the supplied Mapping "smap", so modify the other
   supplied arguments so that they refer to the single output of "map". */
         astInvert( map );
         lbnd_in += ins[ 0 ];
         ubnd_in += ins[ 0 ];
         coord_out = 0;

/* If the output was fed by more than one input, annul the split mapping
   and use the supplied nmapping. */
      } else {
         (void) astAnnul( map );
         map = astClone( smap );
      }
      ins = astFree( ins );

/* If the supplied Mapping could not be split, use the supplied nmapping. */
   } else {
      map = astClone( smap );
   }

/* Check the Mapping is defined in both directions. */
   if( astGetTranForward( map ) && astGetTranInverse( map ) ) {

/* Allocate resources. */
      nin = astGetNin( map );
      nout = astGetNout( map );
      xl = astMalloc( sizeof( double )*(size_t) nin );
      pset1 = astPointSet( NP, nin, "", status );
      ptr1 = astGetPoints( pset1 );
      pset2 = astPointSet( NP, nout, "", status );
      ptr2 = astGetPoints( pset2 );

/* Call astMapBox in a new error reporting context. */
      boxok = 0;
      if( astOK ) {

/* Temporarily switch off error reporting so that no report is made if
   astMapBox cannot find a bounding box (which can legitimately happen with
   some non-linear Mappings). */
         oldrep = astReporting( 0 );

/* Find the upper and lower bounds on the specified Mapping output. This also
   returns the input coords of a point at which the required output has its
   lowest value. */
         astMapBox( map, lbnd_in, ubnd_in, 1, coord_out, &lbnd_out, &ubnd_out,
                    xl, NULL );

/* If the box could not be found, clear the error status and pass on. */
         if( !astOK ) {
            astClearStatus;

/* If the box was found OK, flag this and check if the bounds are equal.
   If so we cannot use them. In this case create new bounds. */
         } else {
            boxok = 1;
            if( astEQUAL( lbnd_out, ubnd_out ) ) {
               m = 0.5*( lbnd_out + ubnd_out );
               if( fabs( m ) > 1.0E-15 ) {
                  lbnd_out = 0.9*m;
                  ubnd_out = 1.1*m;
               } else {
                  lbnd_out = -1.0;
                  ubnd_out = 1.0;
               }
            }
         }

/* Re-instate error reporting. */