Re: [rsf-devel] Issue with sf_floatwrite outputting complex

[Carlos A. da C. F. <c.d...@gm...>]

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


On Fri, Jun 30, 2023 at 1:05 PM Aaron Girard <aar...@gm...>
wrote:

> 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
> _______________________________________________
> RSF-devel mailing list
> RSF...@li...
> https://lists.sourceforge.net/lists/listinfo/rsf-devel
>