C/C++ + SLLIB + SFITSIOによるデータ解析講習会 [Part2] 解答例

山内@みさと天文台/JAXA

2013年7月 初版
2015年11月 第2版
会場: 天文データセンター


課題1

  1. #include <stdio.h>
    #include <math.h>
    
    namespace safe    /* ネームスペース「safe」の開始 */
    {
    
    inline static double acos( double x )
    {
        const double d = 0.0000001;    /* 許容する誤差 */
        if ( x < -1 ) {
            if ( x < -1 - d ) ::fprintf(stderr,"[WARNING] %s: arg is %.15g\n",
                                        __PRETTY_FUNCTION__, x);
            x = -1;
        }
        else if ( 1 < x ) {
            if ( 1 + d < x ) ::fprintf(stderr,"[WARNING] %s: arg is %.15g\n",
                                       __PRETTY_FUNCTION__, x);
            x = 1;
        }
        return ::acos(x);      /* 「::」を書き忘れると自身を呼び出してしまうので注意 */
    }
    
    inline static float acos( float x )
    {
        const double d = 0.0000001;    /* 許容する誤差 */
        if ( x < -1 ) {
            if ( x < -1 - d ) ::fprintf(stderr,"[WARNING] %s: arg is %.15g\n",
                                        __PRETTY_FUNCTION__, x);
            x = -1;
        }
        else if ( 1 < x ) {
            if ( 1 + d < x ) ::fprintf(stderr,"[WARNING] %s: arg is %.15g\n",
                                       __PRETTY_FUNCTION__, x);
            x = 1;
        }
        return ::acosf(x);
    }
    
    }                 /* ネームスペース「safe」の終了 */
    
    int main()
    {
        ::printf("acos = %g\n", ::acos(1.1));
        ::printf("safe acos = %g\n", safe::acos((double)1.1));
        ::printf("safe acos = %g\n", safe::acos((float)1.1));
    
        return 0;
    }
    

課題2

  1. #include <stdio.h>
    #include <math.h>
    
    template <class datatype>
    inline static double get_mean( const datatype vals[], size_t n )
    {
        double sum = 0;
        size_t i;
        for ( i=0 ; i < n ; i++ ) sum += vals[i];
        return sum / (double)n;
    }
    
    int main()
    {
        double v0[] = {10,20,30,40,50,60,70,80,90,100};
        const size_t n_v0 = sizeof(v0) / sizeof(*(v0));
        float v1[] = {1,2,3,4,5,6,7,8,9,10};
        const size_t n_v1 = sizeof(v0) / sizeof(*(v0));
    
        printf("mean of v0 = %g\n", get_mean(v0,n_v0));
        printf("mean of v1 = %g\n", get_mean(v1,n_v1));
    
        return 0;
    }
    

課題3

  1. 2. をご覧ください.
  2. #include <sli/stdstreamio.h>
    #include <sli/digeststreamio.h>
    #include <sli/tarray_tstring.h>
    #include <sli/tstring.h>
    #include <math.h>
    using namespace sli;
    
    int main( int argc, char *argv[] )
    {
        int ret_status = -1;
        stdstreamio sio;                            /* stdin, stdout and stderr */
        digeststreamio f_in;                        /* for local file */
    
        if ( 1 < argc ) {
    
            const char *filename = argv[1];
            tstring line;                           /* buffer for a line */
            tarray_tstring elms;
            
            /* open a stream */
            if ( f_in.open("r", filename) < 0 ) {
                sio.eprintf("[ERROR] cannot open: %s\n",filename);
                goto quit;
            }
    
            /* read text line by line */
            while ( (line=f_in.getline()) != NULL ) {
                double ra, dec;
    
                line.chomp();                               /* erase "\n" */
    
                elms.split(line, " ", false);
    
                ra = M_PI * elms[1].atof() / 180.0;
                dec = M_PI * elms[2].atof() / 180.0;
    
                if ( elms.length() == 3 ) {
                    sio.printf("%s %s %s %g %g %g\n",
                               elms[0].cstr(),elms[1].cstr(),elms[2].cstr(),
                               cos(ra)*cos(dec), sin(ra)*cos(dec), sin(dec));
                }
            }
    
            /* close stream */
            f_in.close();
    
        }
    
        ret_status = 0;
     quit:
        return ret_status;
    }
    
    実行結果(最初の4行):
    M1 083.63308 +22.01450 0.10281 0.921371 0.374841
    M2 323.36258 -00.82325 0.802345 -0.596687 -0.0143679
    M3 205.54842 +28.37728 -0.793808 -0.379451 0.475275
    M4 245.89675 -26.52575 -0.365393 -0.816723 -0.4466
    

課題4

  1. #include <sli/stdstreamio.h>
    #include <sli/fitscc.h>
    using namespace sli;
    
    int main( int argc, char *argv[] )
    {
        int return_status = -1;
        stdstreamio sio;                                    /* standard I/O */
        fitscc fits;                                        /* FITS object  */
    
        /* reading a FITS file */
        if ( 1 < argc ) {
    
            long i, j, npix;
            double sum;
    
            const char *in_file = argv[1];
            if ( fits.read_stream(in_file) < 0 ) {
                sio.eprintf("[ERROR] fits.read_stream(\"%s\") failed\n",in_file);
                goto quit;
            }
    
            fits_image &prim = fits.image("Primary");
    
            prim.convert_type(FITS::FLOAT_T);
    
            fits::float_t *const *ptr_2d;
            ptr_2d = prim.float_t_ptr_2d(true);
    
            sum = 0;
            npix = 0;
            for ( i=0 ; i < prim.row_length() ; i++ ) {
                for ( j=0 ; j < prim.col_length() ; j++ ) {
                    double v = ptr_2d[i][j];
                    if ( isfinite(v) ) {
                        sum += v;
                        npix ++;
                    }
                }
            }
    
            sio.printf("mean = %g\n", sum/npix);
    
        }
    
        return_status = 0;
     quit:
        return return_status;
    }
    
    実行結果:
    mean = 962.69
    

課題5

  1. #include <sli/stdstreamio.h>
    #include <sli/fitscc.h>
    #include <sli/mdarray_statistics.h>
    using namespace sli;
    
    int main( int argc, char *argv[] )
    {
        int return_status = -1;
        stdstreamio sio;                                    /* standard I/O */
        fitscc fits;                                        /* FITS object  */
    
        /* reading a FITS file */
        if ( 1 < argc ) {
    
            const char *in_file = argv[1];
            if ( fits.read_stream(in_file) < 0 ) {
                sio.eprintf("[ERROR] fits.read_stream(\"%s\") failed\n",in_file);
                goto quit;
            }
    
            fits_image &prim = fits.image("Primary");
    
            prim.convert_type(FITS::FLOAT_T);
    
            mdarray_float &prim_arr = prim.float_array();
    
            /* get mean, variance, skewness, kurtosis */
            double stddev;
            mdarray_double moment = md_moment(prim_arr, 
                                              false, NULL, &stddev);
    
            sio.printf("ALL:\n");
            sio.printf("mean = %g\n", moment[0]);
            sio.printf("stddev = %g\n", stddev);
    
            if ( 0 < prim.header_value_length("BIASSEC") ) {
    
                const char *biassec = prim.header("BIASSEC").svalue();
                
                /* get mean, variance, skewness, kurtosis */
                mdarray_double moment = md_moment(prim_arr.sectionf(biassec), 
                                                  false, NULL, &stddev);
    
                sio.printf("OVERSCAN:\n");
                sio.printf("mean = %g\n", moment[0]);
                sio.printf("stddev = %g\n", stddev);
    
            }
    
        }
    
        return_status = 0;
     quit:
        return return_status;
    }
    
    実行結果:
    ALL:
    mean = 962.69
    stddev = 1348.48
    OVERSCAN:
    mean = 920.816
    stddev = 5.93222