[rsf-devel] Issue with sf_floatwrite outputting complex

[Aaron Girard]

Hello,

I am working on a piece of code and am running into a problem concerning float versus complex output. The input is a complex array and the output should be a float array of the magnitude, but it is not.

I wrote a test that takes an input in complex values, calculates the complex magnitude, and should output a float array of the same size. It compiles, but when run the output file is:

esize=8 type=complex 

instead of:

esize=4 type=float, 

and the file says it is 50% of it's expected length.

Below is the example code, and an example of running it. Please let me know if yu have any suggestions that could help me solve this problem.

/* test input complex output a float
  */

#include <math.h>
#include <float.h>
#include <complex.h>

#include <rsf.h>

int main(int argc, char* argv[])
{
    int ny,iy;
    int nx,ix;

    float dy,oy;
    float dx,ox;

    sf_complex **data;
    float **disp;

    sf_axis ax,ay;
    sf_file in,out;

    sf_init(argc,argv);
    /*------------------------------------------------------------*/
    in = sf_input("in");
    out = sf_output("out");

    /* read input file parameters */
    if (SF_COMPLEX != sf_gettype(in)) sf_error("Need complex input");

    if (!sf_histint(in,"n1",&ny)) sf_error("No n1= in input");
    if (!sf_histfloat(in,"d1",&dy)) sf_error("No d1= in input");
    if (!sf_histfloat(in,"o1",&oy)) oy=0.;

    if (!sf_histint(in,"n2",&nx)) sf_error("No n2= in input");
    if (!sf_histfloat(in,"d2",&dx)) sf_error("No d2= in input");
    if (!sf_histfloat(in,"o2",&ox)) ox=0.;

    /* output file parameters */
    ax = sf_maxa(nx,ox,dx);
    sf_oaxa(out,ax,1);

    ay = sf_maxa(ny,oy,dy);
    sf_oaxa(out,ay,2);

    /* memory allocations */
    data     = sf_complexalloc2(nx,ny);
    disp     = sf_floatalloc2(nx,ny);


    for (iy = 0; iy < ny; iy++) {
        for (ix = 0; ix < nx; ix++) {
            data[iy][ix]     = 0.;

       }
    }
    sf_complexread(data[0],ny*nx,in);

    for (iy = 0; iy < ny; iy++) {
        for (ix = 0; ix < nx; ix++) {
            disp[iy][ix] = 0.;
        }
    }
    for (iy = 0; iy < ny; iy++) {
        for (ix = 0; ix < nx; ix++) {
            disp[iy][ix] = sqrt(crealf(data[iy][ix])*crealf(data[iy][ix]) + cimagf(data[iy][ix])*cimagf(data[iy][ix]));
        }
    }
    /* output dispersion */
    sf_floatwrite(disp[0],nx*ny,out);

    exit(0);
}

 ---------------------------------------------------------------

[agirard@c038 Aaron]$ sfspike n1=100 d1=0.002 k1=15 | sfbandpass flo=1 fhi=90 | sfspray axis=2 n=100 d=5 | sffft1 | $RSFSRC/user/agirard/sftest_comp > foo.rsf
sffft1: using 1 of 1 threads

[agirard@c038 Aaron]$ sfin foo.rsf
foo.rsf:
    in="/projects/cwpshell/scratch/foo.rsf@"
    esize=8 type=complex form=native
    n1=100         d1=5           o1=0          label1="" unit1=""
    n2=51          d2=5           o2=0          label2="" unit2=""
    n3=1           d3=?           o3=?          
5100 elements 40800 bytes
sfin: Actually 20400 bytes, 50% of expected.

___________________

Aaron Girard, PhD

Research Associate 

Department of Geophysics 

Colorado School of Mines

[Carlos Alberto da Costa Filho]

Hi Aaron,

As I understand it, if you don't explicitly set the type of the output (or any other parameter), it will inherit from the input. In this case, a simple

sf_settype(out, SF_FLOAT);

should suffice.

Best,

Carlos

--

Carlos Alberto da Costa Filho, Ph.D.

Senior Research Scientist

_______________________________________________
RSF-devel mailing list
RSF-...@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/rsf-devel