Merge commit '7fde2daeed593684120d75de07598154f3ddaf2c' as 'ast'

1 files changed, 864 insertions, 0 deletions

diff --git a/ast/devtools/simtest.c b/ast/devtools/simtest.c
new file mode 100644
index 0000000..e36634c
--- /dev/null
+++ b/ast/devtools/simtest.c
@@ -0,0 +1,864 @@
+#define astCLASS simtest
+
+#include <stdio.h>
+#include <string.h>
+#include "globals.h"
+#include "object.h"
+#include "fitschan.h"
+#include "mapping.h"
+#include "cmpmap.h"
+#include "error.h"
+#include "winmap.h"
+#include "tranmap.h"
+#include "frameset.h"
+#include "pointset.h"
+#include "f77.h"
+#include "memory.h"
+#include "switchmap.h"
+#include "unitmap.h"
+
+void DisplayList( AstMapping *, int, int, int * );
+void sink( const char * );
+const char *source( void );
+void tran( AstMapping *, int * );
+void diff( AstMapping *, int * );
+void tranm( AstMapping **, int *, int, int, int * );
+
+extern F77_SUBROUTINE(kpg1_asreg)( INTEGER(status) );
+
+
+FILE *fd=NULL;
+
+main(){
+   char file[ 50 ], text[ 100 ];
+   const char *class, *enc0;
+   char *attrs, *kvalue;
+   AstWinMap *winmap;
+   AstFitsChan *fc, *fc2;
+   AstChannel *ch, *channel;
+   AstMapping *map, *map2, *tmap;
+   AstObject *obj, *obj2;
+   int opt, nin, nout, i;
+   DECLARE_INTEGER(lstatus);
+   int *status;
+
+   lstatus = 0;
+   status = &lstatus;
+
+/*   F77_CALL(kpg1_asreg)( INTEGER_ARG(&lstatus) ); */
+
+   printf( "Input FITS or AST file [| fitschan attrs] > ");
+   scanf( "%s", file );
+
+   attrs = strchr( file, '|' );
+   if( attrs ) {
+      *attrs = 0;
+      attrs++;
+   } else {
+      attrs = "Carlin=0";
+   }
+
+   fd = fopen( file, "r" );
+   if( fd ){
+      ch = (AstChannel *) astFitsChan( NULL, sink, attrs, status );
+      if( astOK ) {
+         fc = (AstFitsChan *) ch;
+         while( astOK && fgets( text, 100, fd ) ){
+            astPutFits( fc, text, 0 );
+         }
+         astClearCard( fc );
+
+         enc0 = astGetC( fc, "encoding" );
+         printf("Using %s encoding\n", enc0 );
+         obj = astRead( fc );
+      }
+      if( !obj || !astOK ){
+         printf("Attempting to read file %s as an AST file...\n", file );
+         astClearStatus;
+         ch = astChannel( source, NULL, "", status );
+
+         fclose( fd );
+         fd = NULL;
+         fd = fopen( file, "r" );
+         obj = astRead( ch );
+         if( obj != NULL ) {
+            fc = NULL;
+         } else {
+            astClearStatus;
+            printf("Reverting to FitsChan and using a dummy UnitMap. \n", file );
+            obj = (AstObject *) astUnitMap( 1, "", status );
+            astAnnul( ch );
+            ch = (AstChannel *) fc;
+         }
+
+      }
+      fclose( fd );
+      fd = NULL;
+   } else {
+      astError( 1, "Cannot open file '%s'.", status, file );
+   }
+
+   class = obj ? astGetClass( obj ) : NULL;
+   if( obj && class ){
+      printf( "%s %s has been read from file '%s'.\n",
+              strspn( class, "AEIOUaeiou" )?"An":"A", class, file );
+   }
+
+   if( obj && astIsAFitsChan( obj ) ) {
+      fc = (AstFitsChan *) obj;
+      astClearCard( fc );
+      enc0 = astGetC( fc, "encoding" );
+      printf("Using %s encoding\n", enc0 );
+      obj = astRead( fc );
+   }
+
+   if( obj && astIsAMapping( obj ) ){
+      map = (AstMapping *) obj;
+
+      while( astOK ){
+         printf("\n"
+                "0 - Quit\n"
+                "1 - Transform points\n"
+                "2 - Simplify\n"
+                "3 - Display\n"
+                "4 - Write\n"
+                "5 - Shift and scale\n"
+                "6 - Set mapping attributes\n"
+                "7 - Get mapping attributes\n"
+                "8 - Dump the FitsChan\n"
+                "9 - Differentiate\n"
+                "10 - Find a FITS keyword\n"
+                "\n"
+                "Option > " );
+         opt = 0;
+         scanf( "%d", &opt );
+
+         if( opt == 0 ){
+            break;
+
+         } else if( opt == 1 ){
+            tran( map, status );
+
+         } else if( opt == 2 ){
+            map2 = astSimplify( map );
+            map = astAnnul( map );
+            map = map2;
+            obj = (AstObject *) map2;
+
+         } else if( opt == 3 ){
+            if( astIsAFrameSet( map ) ) {
+               tmap = astGetMapping( map, AST__BASE, AST__CURRENT );
+            } else {
+               tmap = astClone( map );
+            }
+            DisplayList( tmap, 1, astGetInvert( tmap ), status );
+            tmap = astAnnul( tmap );
+
+         } else if( opt == 4 ){
+            printf( "\n  Output FITS file > ");
+            scanf( "%s", file );
+            fd = fopen( file, "w" );
+            if( fd ){
+               if( fc ) {
+                  printf("  Encoding (orig=%s)> ",enc0 );
+                  scanf( "%s", text );
+                  if( !strcmp( enc0, text ) && fc != NULL ){
+                     printf("Using existing FitsChan.\n");
+                     fc2 = astClone( fc );
+                     fc = astAnnul( fc );
+                     ch = NULL;
+                     astClearCard( fc2 );
+                     astSet( fc2, "encoding=%s", status, text );
+                  } else {
+                     printf("Creating new FitsChan.\n");
+                     fc2 = astFitsChan( NULL, sink, "encoding=%s", status, text );
+                     astPutFits( fc2, "NAXIS1  = 1", 1 );
+                     astPutFits( fc2, "NAXIS2  = 1", 1 );
+                     astPutFits( fc2, "NAXIS3  = 1", 1 );
+                     astPutFits( fc2, "NAXIS4  = 1", 1 );
+                  }
+               } else {
+                  printf("  Encoding for new FitsChan > " );
+                  scanf( "%s", text );
+                  fc2 = astFitsChan( NULL, sink, "encoding=%s", status, text );
+                  astPutFits( fc2, "NAXIS1  = 1", 1 );
+                  astPutFits( fc2, "NAXIS2  = 1", 1 );
+                  astPutFits( fc2, "NAXIS3  = 1", 1 );
+                  astPutFits( fc2, "NAXIS4  = 1", 1 );
+               }
+               if( !astWrite( fc2, obj ) ){
+                  printf("No object written to file '%s'.\n", file );
+               }
+               fc2 = astAnnul( fc2 );
+               fclose( fd );
+               fd = NULL;
+            } else {
+               astError( 1, "Cannot open file '%s'.", status, file );
+            }
+
+         } else if( opt == 5 ){
+            nin = astGetNin( map );
+            winmap = astWinMap( nin, NULL, NULL, NULL, NULL, "", status );
+
+            printf( "\n   Enter %d axis shift values > ", nin );
+            for( i = 0; i < nin; i++ ) scanf( "%lg", (winmap->a) + i );
+
+            printf( "\n   Enter %d axis scale values > ", nin );
+            for( i = 0; i < nin; i++ ) scanf( "%lg", (winmap->b) + i );
+
+            if( !strcmp( class, "FrameSet" ) ){
+               astInvert( winmap );
+               astRemapFrame( obj, AST__BASE, winmap );
+            } else {
+               map2 = (AstMapping *) astCmpMap( winmap, map, 1, "", status );
+               map = astAnnul( map );
+               map = map2;
+               obj = (AstObject *) map2;
+            }
+            winmap = astAnnul( winmap );
+
+         } else if( opt == 6 ) {
+            printf("Enter attribute specifications: ");
+            scanf( "%s", text );
+            astSet( map, text, status );
+            if( !astOK ) astClearStatus;
+
+         } else if( opt == 7 ) {
+            printf("Enter attribute name: ");
+            scanf( "%s", text );
+            printf( "%s = %s\n", text, astGetC( map, text ) );
+            if( !astOK ) astClearStatus;
+
+         } else if( opt == 8 ) {
+            printf( "\n  Output text file > ");
+            scanf( "%s", file );
+            fd = fopen( file, "w" );
+            if( fd ){
+               channel = astChannel( NULL, sink, "Full=1", status );
+               astDump( ch, channel );
+               channel = astAnnul( channel );
+               fclose( fd );
+               fd = NULL;
+            } else {
+               printf("Failed to open output file '%s'.\n", file );
+               perror(" ");
+            }
+
+         } else if( opt == 9 ){
+            diff( map, status );
+
+         } else if( opt == 10 ){
+            if( fc ) {
+               printf("Enter FITS keyword: ");
+               scanf( "%s", text );
+               if( astGetFitsS( fc, text, &kvalue ) ) {
+                  printf("%s = '%s'\n", text, kvalue );
+               } else {
+                  printf("%s not found\n", text );
+               }
+            } else {
+               printf("No FITS keywords are available\n", text );
+            }
+         }
+      }
+
+   } else {
+      astError( 1, "Object should be a Mapping).", status );
+   }
+
+   if( ch ) ch = astAnnul( ch );
+   obj = astAnnul( obj );
+
+}
+
+void DisplayList( AstMapping *map0, int series, int invert, int *status ){
+   int nmap, i, inv, res, junk;
+   AstMapping *mi, *mf, **map_list, *m, *m2, *map, *smap, *tm1, *tm2;
+   AstCmpMap *cm;
+   AstTranMap *tm;
+   AstChannel *channel;
+   int invi, invf, *invert_list, opt, old_inv, nin, in[20], *out, nout;
+   int isswitch, tfwd, tinv;
+   const char *class;
+   char file[ 50 ];
+   double dv;
+
+   nmap = 0;
+   map_list = NULL;
+   invert_list = NULL;
+   map = astClone( map0 );
+
+   isswitch = astIsASwitchMap( map );
+   if( isswitch ) {
+      astSwitchList( (AstSwitchMap *) map, invert, &nmap, &map_list,
+                     &invert_list, status );
+   } else {
+      astMapList( map, series, invert, &nmap, &map_list, &invert_list );
+   }
+
+   while( 1 ){
+      if( isswitch ) {
+         printf("\nSwitchMap: (fwd selector, inv selector, routes )\n" );
+         printf("  -3 - <transform points>\n");
+      } else {
+         printf("\nMappings in %s:\n", series?"series":"parallel" );
+         printf("  -8 - <invert>\n");
+         printf("  -7 - <write list to text file>\n");
+         printf("  -6 - <split>\n");
+         printf("  -5 - <simplify whole list>\n");
+         printf("  -4 - <merge a single mapping>\n");
+         printf("  -3 - <transform points>\n");
+         printf("  -2 - <show to text file>\n");
+         printf("  -1 - <show to screen>\n");
+      }
+      printf("  0 - <up>\n");
+      for( i = 0; i < nmap; i++ ){
+         m = map_list[ i ];
+         class = m ? astGetClass( m ) :"(NULL)";
+         printf( "  %d - %s\n", i + 1, class );
+      }
+      printf("\n  Option > ", 0 );
+      opt = 0;
+      scanf("%d", &opt );
+
+      if ( opt == -8 ) {
+
+         for( i = 0; i < nmap; i++ ) map_list[ i ] = astAnnul( map_list[ i ] );
+         map_list = astFree( map_list );
+         invert_list = astFree( invert_list );
+         nmap = 0;
+
+         invert = ( invert == 0 );
+         if( isswitch ) {
+            astSwitchList( (AstSwitchMap *) map, invert, &nmap,
+                            &map_list, &invert_list, status );
+         } else {
+            astMapList( map, series, invert, &nmap, &map_list, &invert_list );
+         }
+
+      } else if ( opt == -7 && !isswitch ) {
+         printf( "\n  Output text file > ");
+         scanf( "%s", file );
+         fd = fopen( file, "w" );
+         if( fd ){
+            channel = astChannel( NULL, sink, "", status );
+            for( i = 0; i < nmap; i++ ) {
+               inv = astGetInvert( map_list[ i ] );
+               astSetInvert(  map_list[ i ], invert_list[ i ] );
+               astDump( map_list[ i ], channel );
+               astSetInvert(  map_list[ i ], inv );
+            }
+            channel = astAnnul( channel );
+            fclose( fd );
+            fd = NULL;
+         } else {
+            printf("Failed to open output file.\n" );
+         }
+
+      } else if( opt == -6 && !isswitch ) {
+         old_inv = astGetInvert( map );
+         astSetInvert( map, invert );
+
+         printf( "\n   Pick how many input axes (max=%d)? > ",
+                 astGetNin( map ) );
+         scanf( "%lg", &dv );
+         nin = dv;
+
+         printf( "\n   Enter %d zero-based axis indices > ", nin );
+         for( i = 0; i < nin; i++ ) {
+            scanf( "%lg", &dv );
+            in[ i ] = dv;
+         }
+
+         out = astMapSplit( map, nin, in, &smap );
+         astSetInvert( map, old_inv );
+
+         if( out ) {
+            printf( "Corresponding outputs are: " );
+            nout = astGetNout( smap );
+            for( i = 0; i < nout; i++ ) printf( "%d ", out[ i ] );
+            printf("\n");
+            out = astFree( out );
+            map = astAnnul( map );
+            map = smap;
+         } else {
+            printf( "Selected inputs do not correspond to a subset of outputs\n");
+         }
+
+         for( i = 0; i < nmap; i++ ) map_list[ i ] = astAnnul( map_list[ i ] );
+         map_list = astFree( map_list );
+         invert_list = astFree( invert_list );
+
+         nmap = 0;
+         astMapList( map, series, invert, &nmap, &map_list, &invert_list );
+
+      } else if( opt == -5 && !isswitch ) {
+         old_inv = astGetInvert( map );
+         astSetInvert( map, invert );
+         smap = astSimplify( map );
+         if( smap == map ) printf("\n  No simplification performed\n" );
+         astSetInvert( map, old_inv );
+         map = astAnnul( map );
+         map = smap;
+
+         for( i = 0; i < nmap; i++ ) map_list[ i ] = astAnnul( map_list[ i ] );
+         map_list = astFree( map_list );
+         invert_list = astFree( invert_list );
+
+         nmap = 0;
+         astMapList( map, series, invert, &nmap, &map_list, &invert_list );
+
+      } else if( opt == -4 && !isswitch ) {
+         printf("\n  Mapping index > ", 0 );
+         scanf("%d", &i );
+         i--;
+         res = astMapMerge( map_list[ i ], i, series, &nmap,
+                            &map_list, &invert_list );
+         if( res != -1 ) {
+            printf( "Index of first modified Mapping: %d\n", res + 1 );
+         } else {
+            printf( "No mappings modified\n" );
+         }
+
+      } else if( opt == -3 && series && !isswitch ) {
+
+         if( invert == astGetInvert( map ) ) {
+            tfwd = astGetTranForward( map );
+            tinv = astGetTranInverse( map );
+         } else {
+            tinv = astGetTranForward( map );
+            tfwd = astGetTranInverse( map );
+         }
+
+         if( !tfwd ) {
+            printf( "\n   No forward transformation available\n" );
+
+         } else if( !tinv ) {
+            printf( "\n   No inverse transformation available\n" );
+            tranm( map_list, invert_list, nmap, 0, status );
+
+         } else {
+            tranm( map_list, invert_list, nmap, 1, status );
+
+         }
+
+      } else if( opt == -3 && isswitch ) {
+
+         if( invert == astGetInvert( map ) ) {
+            tfwd = astGetTranForward( map );
+            tinv = astGetTranInverse( map );
+         } else {
+            tinv = astGetTranForward( map );
+            tfwd = astGetTranInverse( map );
+         }
+
+         if( !tfwd ) {
+            printf( "\n   No forward transformation available\n" );
+
+         } else if( !tinv ) {
+            printf( "\n   No inverse transformation available\n" );
+            tranm( &map, &invert, 1, 0, status );
+
+         } else {
+            tranm( &map, &invert, 1, 1, status );
+
+         }
+
+      } else if( opt == -3 && !series && !isswitch ) {
+         old_inv = astGetInvert( map );
+         astSetInvert( map, invert );
+         tran( map, status );
+         astSetInvert( map, old_inv );
+
+      } else if ( opt == -2 && !isswitch ) {
+
+         for( i = 0; i < nmap; i++ ) {
+            inv = astGetInvert( map_list[ i ] );
+            astSetInvert(  map_list[ i ], invert_list[ i ] );
+            if( i == 0 ) {
+               tm1 = astClone( map_list[ 0 ] );
+            } else {
+               tm2 = (AstMapping *) astCmpMap( tm1, map_list[ i ], series, "", status );
+               tm1 = astAnnul( tm1 );
+               tm1 = tm2;
+            }
+            astSetInvert(  map_list[ i ], inv );
+         }
+
+         printf( "\n  Output text file > ");
+         scanf( "%s", file );
+         fd = fopen( file, "w" );
+         if( fd ){
+            channel = astChannel( NULL, sink, "", status );
+            old_inv = astGetInvert( tm1 );
+            astSetInvert( tm1, invert );
+            astDump( tm1, channel );
+            astSetInvert( tm1, old_inv );
+            channel = astAnnul( channel );
+            fclose( fd );
+            fd = NULL;
+         } else {
+            printf("Failed to open output file.\n" );
+         }
+         tm1= astAnnul( tm1 );
+
+      } else if ( opt == -1 && !isswitch ) {
+
+
+         for( i = 0; i < nmap; i++ ) {
+            inv = astGetInvert( map_list[ i ] );
+            astSetInvert(  map_list[ i ], invert_list[ i ] );
+            if( i == 0 ) {
+               tm1 = astClone( map_list[ 0 ] );
+            } else {
+               tm2 = (AstMapping *) astCmpMap( tm1, map_list[ i ],
+                                               series, "", status );
+               tm1 = astAnnul( tm1 );
+               tm1 = tm2;
+            }
+            astSetInvert(  map_list[ i ], inv );
+         }
+
+         channel = astChannel( NULL, NULL, "Full=1", status );
+         old_inv = astGetInvert( tm1);
+         astSetInvert( tm1, invert );
+         astDump( tm1, channel );
+         astSetInvert( tm1, old_inv );
+         channel = astAnnul( channel );
+         tm1= astAnnul( tm1 );
+
+      } else if( opt > 0 && opt <= nmap && map_list[ opt - 1 ] ){
+         m2 = map_list[ opt - 1 ];
+         inv = invert_list[ opt - 1 ];
+         class = astGetClass( m2 );
+
+         if( m2 && astIsAFrameSet( m2 ) ){
+            m = astGetMapping( m2, AST__BASE, AST__CURRENT );
+            class = astGetClass( m );
+
+         } else if( astIsATranMap( m2 ) ){
+
+            astDecompose( m2, &mf, &mi, &invf, &invi, &junk );
+
+            if( inv ) {
+               m = mi;
+               inv = invi;
+            } else {
+               m = mf;
+               inv = invf;
+            }
+            class = astGetClass( m );
+
+         } else {
+            m = astClone( m2 );
+         }
+
+         if( strcmp( class, "CmpMap" ) ){
+            DisplayList( m, 1, inv, status );
+         } else {
+            cm = (AstCmpMap *) m;
+            DisplayList( m, cm->series, inv, status );
+         }
+      } else {
+         break;
+      }
+   }
+
+   for( i = 0; i < nmap; i++ ) map_list[ i ] = astAnnul( map_list[ i ] );
+   map_list = astFree( map_list );
+   invert_list = astFree( invert_list );
+   nmap = 0;
+
+}
+
+
+void sink( const char *text ){
+   if( fd && text ) fprintf( fd, "%s\n", text );
+}
+
+const char *source( void ){
+   static char text[2000];
+   static char *ptext;
+   if( fd ) {
+      ptext = fgets( text, 2000, fd );
+   } else {
+      ptext = NULL;
+   }
+   return ptext;
+}
+
+void tran( AstMapping *map, int *status ){
+   char text[2000];
+   AstFrame *infrm, *outfrm;
+   AstPointSet *pset1, *pset2;
+   double **ptr1, **ptr2;
+   int nin, nout, i, ok;
+
+   if( !astGetTranForward( map ) ) {
+      printf( "\n   No forward transformation available\n" );
+      return;
+   }
+
+   if( !astGetTranInverse( map ) ) {
+      printf( "\n   No inverse transformation available\n" );
+   }
+
+   nin = astGetNin( map );
+   nout = astGetNout( map );
+   pset1 = astPointSet( 1, nin, "", status );
+   ptr1 = astGetPoints( pset1 );
+   pset2 = astPointSet( 1, nout, "", status );
+   ptr2 = astGetPoints( pset2 );
+
+   if( astIsAFrameSet( map ) ) {
+      infrm = astGetFrame( map, AST__BASE );
+      outfrm = astGetFrame( map, AST__CURRENT );
+   } else if( astIsAFrame( map ) ) {
+      infrm = (AstFrame *) map;
+      outfrm = (AstFrame *) map;
+   } else {
+      infrm = NULL;
+      outfrm = NULL;
+   }
+
+   if( astOK ){
+      ok = 0;
+      if( infrm != NULL ) {
+         ok = 1;
+         for( i = 0; i < nin; i++ ) {
+            printf( "\n   Enter %s value (or <bad>) > ", astGetLabel( infrm, i ) );
+            scanf( "%s", text );
+            if( !astUnformat( infrm, i, text, ptr1[ i ] ) ) {
+               ok = 0;
+               printf( "Bad formatted value supplied! Now try unformatted axis values...\n");
+               break;
+            }
+         }
+      }
+
+      if( !ok ) {
+         printf( "\n   Enter %d axis values (BAD=999) > ", nin );
+         for( i = 0; i < nin; i++ ) {
+            scanf( "%lg", ptr1[ i ] );
+            if( ptr1[ i ][ 0 ] == 999.0 ) ptr1[ i ][ 0 ] = AST__BAD;
+         }
+      }
+
+      astTransform( map, pset1, 1, pset2 );
+      printf( "\n   Forward: [ " );
+      for( i = 0; i < nin; i++ ) {
+         if( ptr1[ i ][ 0 ] == AST__BAD ) {
+            printf( "999.0 " );
+         } else {
+            printf( "%.*g ", DBL_DIG, ptr1[ i ][ 0 ] );
+         }
+      }
+
+      printf( "] -> [ " );
+      for( i = 0; i < nout; i++ ) {
+         if( ptr2[ i ][ 0 ] == AST__BAD ) {
+            printf( "999.0 " );
+         } else {
+            printf( "%.*g ", DBL_DIG, ptr2[ i ][ 0 ] );
+         }
+      }
+      printf( "]\n" );
+
+      if( infrm != NULL && outfrm != NULL ) {
+         printf( "            ( " );
+         for( i = 0; i < nin; i++ ) {
+            printf( "%s ", astFormat( infrm, i, ptr1[ i ][ 0 ] ) );
+         }
+         printf( ") -> ( " );
+         for( i = 0; i < nout; i++ ) {
+            printf( "%s ", astFormat( outfrm, i, ptr2[ i ][ 0 ] ) );
+         }
+         printf( ") formatted\n" );
+      }
+
+
+      if( astGetTranInverse( map ) ) {
+         astTransform( map, pset2, 0, pset1 );
+         printf( "\n   Inverse: [ " );
+         for( i = 0; i < nout; i++ ) {
+            if( ptr2[ i ][ 0 ] == AST__BAD ) {
+               printf( "999.0 " );
+            } else {
+               printf( "%.*g ", DBL_DIG, ptr2[ i ][ 0 ] );
+            }
+         }
+
+         printf( "] -> [ " );
+         for( i = 0; i < nin; i++ ) {
+            if( ptr1[ i ][ 0 ] == AST__BAD ) {
+               printf( "999.0 " );
+            } else {
+               printf( "%.*g ", DBL_DIG, ptr1[ i ][ 0 ] );
+            }
+         }
+         printf( "]\n" );
+
+
+         if( infrm != NULL && outfrm != NULL ) {
+            printf( "            ( " );
+            for( i = 0; i < nout; i++ ) {
+               printf( "%s ", astFormat( outfrm, i, ptr2[ i ][ 0 ] ) );
+            }
+            printf( ") -> ( " );
+            for( i = 0; i < nin; i++ ) {
+               printf( "%s ", astFormat( infrm, i, ptr1[ i ][ 0 ] ) );
+            }
+            printf( ") formatted\n" );
+         }
+      }
+   }
+   pset1 = astAnnul( pset1 );
+   pset2 = astAnnul( pset2 );
+}
+
+void tranm( AstMapping **map, int *inv_list, int nmap, int inv, int *status ){
+   AstPointSet *pset1, *pset2;
+   double **ptr1, **ptr2;
+   int nin, nout, i, j, inv0;
+
+   inv0 = astGetInvert( map[ 0 ] );
+   astSetInvert( map[ 0 ], inv_list[ 0 ] );
+   nin = astGetNin( map[ 0 ] );
+   astSetInvert( map[ 0 ], inv0 );
+
+   pset1 = astPointSet( 1, nin, "", status );
+   ptr1 = astGetPoints( pset1 );
+
+   if( astOK ){
+      printf( "\n   Enter %d axis values (BAD=999) > ", nin );
+      for( i = 0; i < nin; i++ ) {
+         scanf( "%lg", ptr1[ i ] );
+         if( ptr1[ i ][ 0 ] == 999.0 ) ptr1[ i ][ 0 ] = AST__BAD;
+      }
+
+      printf("Forward...\n");
+      for( j=0; j<nmap; j++ ) {
+         inv0 = astGetInvert( map[ j ] );
+         astSetInvert( map[ j ], inv_list[ j ] );
+
+         nout = astGetNout( map[ j ] );
+         pset2 = astPointSet( 1, nout, "", status );
+         ptr2 = astGetPoints( pset2 );
+
+         astTransform( map[ j ], pset1, 1, pset2 );
+
+         printf( " Mapping %d gave [ ", j+1 );
+         for( i = 0; i < nout; i++ ) {
+            if( ptr2[ i ][ 0 ] == AST__BAD ) {
+               printf( "999.0 " );
+            } else {
+               printf( "%.*g ", DBL_DIG, ptr2[ i ][ 0 ] );
+            }
+         }
+
+         printf( "]\n" );
+
+         pset1 = astAnnul( pset1 );
+         astSetInvert( map[ j ], inv0 );
+
+         pset1 = pset2;
+         ptr1 = ptr2;
+         nin = nout;
+      }
+
+      printf("Inverse...\n");
+
+      for( j=nmap - 1; j>=0 && inv; j-- ) {
+         inv0 = astGetInvert( map[ j ] );
+         astSetInvert( map[ j ], inv_list[ j ] );
+
+         nin = astGetNin( map[ j ] );
+         pset2 = astPointSet( 1, nin, "", status );
+         ptr2 = astGetPoints( pset2 );
+
+         astTransform( map[ j ], pset1, 0, pset2 );
+
+         printf( " Mapping %d gave [ ", j+1 );
+         for( i = 0; i < nin; i++ ) {
+            if( ptr2[ i ][ 0 ] == AST__BAD ) {
+               printf( "999.0 " );
+            } else {
+               printf( "%.*g ", DBL_DIG, ptr2[ i ][ 0 ] );
+            }
+         }
+         printf( "]\n" );
+
+         pset1 = astAnnul( pset1 );
+         astSetInvert( map[ j ], inv0 );
+
+         pset1 = pset2;
+         ptr1 = ptr2;
+         nout = nin;
+      }
+
+      pset1 = astAnnul( pset1 );
+
+   }
+}
+
+
+
+
+void diff( AstMapping *map, int *status ){
+   char text[200];
+   AstFrame *infrm, *outfrm;
+   int nin, nout, i, ok, ax1, ax2;
+   double at[2000];
+
+
+   nin = astGetNin( map );
+   nout = astGetNout( map );
+
+   if( astIsAFrameSet( map ) ) {
+      infrm = astGetFrame( map, AST__BASE );
+      outfrm = astGetFrame( map, AST__CURRENT );
+   } else if( astIsAFrame( map ) ) {
+      infrm = (AstFrame *) map;
+      outfrm = (AstFrame *) map;
+   } else {
+      infrm = NULL;
+      outfrm = NULL;
+   }
+
+   if( astOK ){
+      ok = 0;
+      if( infrm != NULL ) {
+         ok = 1;
+         for( i = 0; i < nin; i++ ) {
+            printf( "\n   Enter %s value (or <bad>) > ", astGetLabel( infrm, i ) );
+            scanf( "%s", text );
+            if( !astUnformat( infrm, i, text, at + i ) ) {
+               ok = 0;
+               printf( "Bad formatted value supplied! Now try unformatted axis values...\n");
+               break;
+            }
+         }
+      }
+
+      if( !ok ) {
+         printf( "\n   Enter %d axis values (BAD=999) > ", nin );
+         for( i = 0; i < nin; i++ ) {
+            scanf( "%lg", at + i );
+            if( at[ i ] == 999.0 ) at[ i ] = AST__BAD;
+            printf("(got %.*g)\n", DBL_DIG, at[ i ] );
+         }
+
+      }
+
+      for( i = 0; i < nin; i++ ) {
+         printf("(%d: got %.*g)\n", i, DBL_DIG, at[ i ] );
+      }
+
+      printf( "\n   Enter output to differentiate (zero-based index): > " );
+      scanf( "%d", &ax1 );
+      printf( "\n   Enter input to vary (zero-based index): > " );
+      scanf( "%d", &ax2 );
+      printf( "%.*g \n", DBL_DIG, astRate( map, at, ax1, ax2 ) );
+   }
+}
+