/***/
/* diag.M: diagonal symmetrization.
/***/

/***/
/* Initialization steps:
/*  ARGS: extract symmetrization function.
/***/

if (sliceCode == CODE_ARGS)
   {
    sumMode = argTest( "-sum" );
   };

if (sliceCode < 1)
   {
    exit( 0 ); 
   };

/***/
/* Process each 2D plane:
/***/

if (sliceCode % ySize)
   {
    exit( 0 );
   };

if (xSize != ySize)
   {
    (void) printf( "Symmetrize Error: non-square data.\n" );
    exit( 1 ); 
   };

float symData[xSize*ySize];
 
(void) dReadB( inUnit,  symData, wordLen*xSize*ySize );

for( iy = 0; iy < ySize; iy++ )
   {
    if (!((iy+1) % 16))
       {
        (void) printf( "Slice %.0f of %.0f\n", iy + 1, ySize );
       };

    for( ix = iy + 1; ix < xSize; ix++ )
       {
        loc1 = ix + iy*xSize;
        loc2 = ix*xSize + iy;

        if (sumMode)
           {
            val = 0.5*(symData[loc1] + symData[loc2]);
           }
        else
           {
            if (abs(symData[loc1]) < abs(symData[loc2]))
               {
                val = symData[loc1];
               }
            else
               {
                val = symData[loc2];
               };
           };

        symData[loc1] = val;
        symData[loc2] = val;
       }
   }

(void) dWrite( outUnit, symData, wordLen*xSize*ySize );
