/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:18 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_srft s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#include "p_srft.h"
		/* PARAMETER translations */
#define	ONE	1.e0
#define	PI	3.1415926535897932384e0
#define	ZERO	0.e0
		/* end of PARAMETER translations */
 
 
int main( )
{
	long int j, j1, j2, k, l, ms, n, n2, n4;
	float a[16][16], s[3], sig, sigd, temp;
	static long m[2]={4,4};
	static long nd = 2;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	long *const M = &m[0] - 1;
	float *const S = &s[0] - 1;
		/* end of OFFSET VECTORS */
 
	/*>> 1996-06-19 DRSRFT Krogh  Minor change for C conversion.
	 *>> 1994-10-19 DRSRFT Krogh  Changes to use M77CON
	 *>> 1994-08-09 DRSRFT WVS Remove '0' from format
	 *>> 1993-02-04 DRSRFT CLL
	 *>> 1989-05-07 DRSRFT FTK, CLL
	 *>> 1989-05-04 DRSRFT FTK, CLL
	 *     Demo driver for SRFT -- Multi-dimensional real Fourier transform
	 *     ------------------------------------------------------------------
	 *--S replaces "?": DR?RFT, ?RFT
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------
	 *                         Start of code -- Construct A */
	n = ipow(2,M[1]);
	n2 = n/2;
	n4 = n2/2;
	sigd = PI/n2;
	for (j1 = 1; j1 <= n; j1++)
	{
		for (j2 = 1; j2 <= n; j2++)
		{
			a[j2 - 1][j1 - 1] = ZERO;
			if (labs( j1 - n2 - 1 ) + labs( j2 - n2 - 1 ) <= n4)
				a[j2 - 1][j1 - 1] = ONE;
		}
	}
	/*     ------------------------------------------------------------------
	 *                Compute Fourier transform and apply sigma factors */
	ms = 0;
	srft( (float*)a, 'A', m, nd, &ms, s );
	for (j1 = 1; j1 <= n; j1 += 2)
	{
		a[n2][j1 - 1] = ZERO;
		a[n2][j1] = ZERO;
		for (j2 = 1; j2 <= n2; j2++)
		{
			sig = ONE;
			if (j1 == 1)
			{
				/*                            No change in SIG due to J */
				if (j2 != 1)
				{
					j = j2 - 1;
					k = 1;
				}
				else
				{
					a[0][1] = ZERO;
					goto L_40;
				}
			}
			else
			{
				j = j1/2;
				k = 0;
			}
			/*                         Get nontrivial sigma factors * SIG */
L_30:
			;
			if (j < n4)
			{
				temp = S[j];
			}
			else if (j == n4)
			{
				temp = ONE;
			}
			else
			{
				temp = S[n2 - j];
			}
			sig = sig*temp/(sigd*(float)( j ));
			if (k == 0)
			{
				if (j2 != 1)
				{
					j = j2 - 1;
					k = 1;
					goto L_30;
				}
			}
			else
			{
				/*                                       Apply sigma factors */
				if (j1 == 1)
				{
					a[n - j2 + 1][0] = ZERO;
					a[n - j2 + 1][1] = ZERO;
				}
				else
				{
					a[n - j2 + 1][j1 - 1] *= sig;
					a[n - j2 + 1][j1] *= sig;
				}
			}
			a[j2 - 1][j1 - 1] *= sig;
			a[j2 - 1][j1] *= sig;
L_40:
			;
		}
	}
	srft( (float*)a, 'S', m, nd, &ms, s );
	printf("\n Smoothed A\n");
	for (l = 1; l <= 9; l++)
	{
		for (n = 1; n <= 9; n++)
		{
			printf("%8.4f", a[n - 1][l - 1]);
		}
		printf("\n");
	}
	exit(0);
} /* end of function */