********************************************************************************
* sinc.c -- Sinc function objects
*
* User routines:
*	Pro_SincFunc()		creates a sinc function object.
*	Pro_SincInfo()		returns parameters of sinc function.
*
* Mark Showalter, PDS Ring-Moon Systems Node, March 1998
*******************************************************************************/
#include <stdio.h>
#include <math.h>
#include "profile.h"
#include "fortran.h"

#define PI 3.141592653589793238462643

/********************
 * Type definitions *
 ********************/

typedef struct ZPRO_SINC_STRUCT {
    RL_INT4	flag;
    RL_FLT8	height, step, halfwidth, shoulder, scale, threshold,
		shoulder_scale;
} ZPRO_SINC;

/**********************
 * Symbol definitions *
 **********************/

#define ZPRO_SINC_FLAG	'sinc'

/************************************
 * Prototypes of internal functions *
 ************************************/

static RL_FLT8 ZPro_SincCalc  RL_PROTO((RL_FLT8 x, RL_VOID *params));
static void    ZPro_PrintSinc RL_PROTO((RL_VOID *params));

/*
********************************************************************************
* EXPORTED USER ROUTINES
********************************************************************************
*$ Component_name:
*	Pro_SincFunc (sinc.c)
*$ Abstract:
*	This routine generates a sinc function object.  A sinc function returns
*	sin(pi*x)/(pi*x), scaled by an envelope function that goes smoothly to
*	zero outside its domain.
*$ Keywords:
*	PROFILE, FUNCTION, PSF
*	C, PUBLIC, SUBROUTINE
*$ Declarations:
*	PRO_OBJECT	*Pro_SincFunc(height, step, halfwidth, shoulder)
*	RL_FLT8         height, step, halfwidth, shoulder;
*$ Inputs:
*	height		peak y-value of function (at origin).
*	step		x-step between origin/zeroes of function.
*	halfwidth	x-value above which y-value is zero.
*	shoulder	x-width of truncation, modeled as a quarter-sinusoid.
*$ Outputs:
*	none
*$ Returns:
*	pointer to a new sinc function object, or NULL on non-fatal error.
*$ Detailed_description:
*	This routine generates a sinc function object.  A sinc function returns
*	sin(pi*x)/(pi*x), scaled by an envelope function that goes smoothly to
*	zero outside its domain.  The x and y coordinates are initially unnamed.
*
*	The envelope function used is as follows:
*		1		|x| < halfwidth - shoulder
*		0		|x| > halfwidth
*		(1 + cos((|x|-(halfwidth-shoulder)) * pi/shoulder))/2
*				otherwise
*	This particular form has the advantage of being smooth throughout the
*	domain.
*$ External_references:
*	Profile toolkit
*$ Side_effects:
*	Memory is allocated.
*$ Examples:
*	This line of code creates a sinc function with height 2 and with zeros
*	spaced by unity, with a full width of 10 and no shoulder transition
*	region.
*
*	sinc = Pro_SincFunc(2., 1., 5., 0.);
*$ Error_handling:
*	Profile toolkit error handling is in effect.
*
*	Conditions raised:
*	RL_MEMORY_ERROR		on memory allocation error.
*$ Limitations:
*	none
*$ Author_and_institution:
*	Mark R. Showalter
*	NASA/Ames Research Center
*$ Version_and_date:
*	1.0: March 1998
*$ Change_history:
*	none
*******************************************************************************/

PRO_OBJECT	*Pro_SincFunc(height, step, halfwidth, shoulder)
RL_FLT8         height, step, halfwidth, shoulder;
{
ZPRO_SINC	*params;
PRO_OBJECT	*func;
RL_INT4		paramsize;

    paramsize = sizeof(ZPRO_SINC);
    params = (ZPRO_SINC *) XRL_Malloc(paramsize);
    if (params == NULL) return NULL;

    params->flag      = ZPRO_SINC_FLAG;
    params->height    = height;
    params->step      = step;
    params->halfwidth = halfwidth;
    params->shoulder  = shoulder;

    params->scale          = PI / step;
    params->threshold      = halfwidth - shoulder;

    params->shoulder_scale = 0.;
    if (params->shoulder != 0.) params->shoulder_scale = PI / shoulder;

    func = Pro_SoftFunc(ZPro_SincCalc, -halfwidth, halfwidth,
		        (RL_VOID *) params, paramsize, TRUE, ZPro_PrintSinc);
    XRL_Free((RL_VOID *) params);

    return func;
}

/*
********************************************************************************
*$ Component_name:
*	Pro_SincInfo (sinc.c)
*$ Abstract:
*	This routine returns information about a sinc function object.
*$ Keywords:
*	PROFILE, FUNCTION, PSF
*	C, PUBLIC, SUBROUTINE
*$ Declarations:
*	RL_BOOL		Pro_SincInfo(sinc, height, step, halfwidth, shoulder)
*	PRO_OBJECT	*sinc;
*	RL_FLT8         *height, *step, *halfwidth, *shoulder;
*$ Inputs:
*	sinc		pointer to a sinc function object.
*$ Outputs:
*	*height		peak y-value of function (at origin);
*			not returned if height==NULL.
*	*step		x-step between origin/zeroes of function;
*			not returned if height==NULL.
*	*halfwidth	x-value above which y-value is zero;
*			not returned if height==NULL.
*	*shoulder	x-width of truncation, modeled as a half-sinusoid;
*			not returned if height==NULL.
*$ Returns:
*	TRUE if this is a sinc function; FALSE otherwise.
*$ Detailed_description:
*	This routine returns information about a sinc function object.
*$ External_references:
*	Profile toolkit
*$ Side_effects:
*	none
*$ Examples:
*	This snippet of code of code prints the step size in a sinc function.
*
*	PRO_OBJECT	*sinc;
*	RL_FLT8		step;
*
*	Pro_SincInfo(sinc, NULL, &step, NULL, NULL);
*	printf("Step size = %g\n", step);
*$ Error_handling:
*	Profile toolkit error handling is in effect.
*
*	Conditions raised:
*	PRO_CLASS_ERROR		if this is not a software function object.
*$ Limitations:
*	none
*$ Author_and_institution:
*	Mark R. Showalter
*	NASA/Ames Research Center
*$ Version_and_date:
*	1.0: March 1998
*$ Change_history:
*	none
*******************************************************************************/

RL_BOOL		Pro_SincInfo(sinc, height, step, halfwidth, shoulder)
PRO_OBJECT	*sinc;
RL_FLT8         *height, *step, *halfwidth, *shoulder;
{
ZPRO_SINC	*params;

    params = (ZPRO_SINC *) XPro_SoftParams(sinc);
    if (params == NULL) return FALSE;

    if (params->flag != ZPRO_SINC_FLAG) return FALSE;

    if (height    != NULL) *height    = params->height;
    if (step      != NULL) *step      = params->step;
    if (halfwidth != NULL) *halfwidth = params->halfwidth;
    if (shoulder  != NULL) *shoulder  = params->shoulder;

    return TRUE;
}

/*
********************************************************************************
* INTERNAL ROUTINES
********************************************************************************
* ZPro_SincCalc(x, params)
*
* This internal routine calculates the value of a sinc function.
*
* Inputs:
*	x		location at which to evaluate function.
*	params		pointer to the ZPRO_SINC data structure.
*
* Return:		value of function at x.
*******************************************************************************/

static RL_FLT8	ZPro_SincCalc(x, params)
RL_FLT8		x;
RL_VOID		*params;
{
ZPRO_SINC	*sincpar;
RL_FLT8		arg, value;

    sincpar = (ZPRO_SINC *) params;

    /* Handle mathematical singularity at origin */
    if (x == 0.) return sincpar->height;

    /* Evaluate sinc function */
    if (x < 0.) x = -x;
    arg = x * sincpar->scale;
    value = sincpar->height * sin(arg)/arg;

    /* If we're not in shoulder region, return */
    if (x <= sincpar->threshold) return value;

    /* Otherwise, scale by shoulder function and return */
    return value * 0.5 *
	(1. + cos((x - sincpar->threshold) * sincpar->shoulder_scale));
}

/*
********************************************************************************
* ZPro_PrintSinc(params)
*
* This internal routine prints information about a sinc function.
*
* Inputs:
*	params		pointer to the ZPRO_SINC data structure.
*******************************************************************************/

static void	ZPro_PrintSinc(params)
RL_VOID		*params;
{
ZPRO_SINC	*sincpar;

    sincpar = (ZPRO_SINC *) params;

    /* Print sinc object info... */
    printf("\nSinc software function parameters...\n");

    printf("   height = %#g\n", sincpar->height);
    printf("     step = %#g\n", sincpar->step);
    printf("halfwidth = %#g\n", sincpar->halfwidth);
    printf(" shoulder = %#g\n", sincpar->shoulder);
}

/*
********************************************************************************
* FORTRAN INTERFACE ROUTINES
********************************************************************************
*$ Component_name:
*	FPro_SincFunc (sinc.c)
*$ Abstract:
*	This routine generates a sinc function object.  A sinc function returns
*	sin(pi*x)/(pi*x), scaled by an envelope function that goes smoothly to
*	zero outside its domain.
*$ Keywords:
*	PROFILE, FUNCTION, PSF
*	FORTRAN, PUBLIC, SUBROUTINE
*$ Declarations:
*	integer*4 function FPro_SincFunc(height, step, halfwidth, shoulder)
*	real*8		height, step, halfwidth, shoulder
*$ Inputs:
*	height		peak y-value of function (at origin).
*	step		x-step between origin/zeroes of function.
*	halfwidth	x-value above which y-value is zero.
*	shoulder	x-width of truncation, modeled as a quarter-sinusoid.
*$ Outputs:
*	none
*$ Returns:
*	FORTRAN pointer to a new sinc function object, or 0 on non-fatal error.
*$ Detailed_description:
*	This routine generates a sinc function object.  A sinc function returns
*	sin(pi*x)/(pi*x), scaled by an envelope function that goes smoothly to
*	zero outside its domain.  The x and y coordinates are initially unnamed.
*
*	The envelope function used is as follows:
*		1		|x| < halfwidth - shoulder
*		0		|x| > halfwidth
*		(1 + cos((|x|-(halfwidth-shoulder)) * pi/shoulder))/2
*				otherwise
*	This particular form has the advantage of being smooth throughout the
*	domain.
*$ External_references:
*	Profile toolkit
*$ Side_effects:
*	Memory is allocated.
*$ Examples:
*	This line of code creates a sinc function with height 2 and with zeros
*	spaced by unity, with a full width of 10 and no shoulder transition
*	region.
*
*	sinc = Pro_SincFunc(2.d0, 1.d0, 5.d0, 0.d0)
*$ Error_handling:
*	Profile toolkit error handling is in effect.
*
*	Conditions raised:
*	RL_MEMORY_ERROR		on memory allocation error.
*$ Limitations:
*	none
*$ Author_and_institution:
*	Mark R. Showalter
*	NASA/Ames Research Center
*$ Version_and_date:
*	1.0: March 1998
*$ Change_history:
*	none
*******************************************************************************/

RL_INT4	FORTRAN_NAME(fpro_sincfunc) (height, step, halfwidth, shoulder)
RL_FLT8	*height, *step, *halfwidth, *shoulder;
{
RL_VOID *ptr;
RL_INT4 index;

    ptr = (RL_VOID *) Pro_SincFunc(*height, *step, *halfwidth, *shoulder);
    if (ptr == NULL) return 0;

    index = FORT_AddPointer(ptr);
    if (index == 0) Pro_FreeObject((PRO_OBJECT *) ptr);

    return index;
}

/*
********************************************************************************
*$ Component_name:
*	FPro_SincInfo (sinc.c)
*$ Abstract:
*	This routine returns information about a sinc function object.
*$ Keywords:
*	PROFILE, FUNCTION, PSF
*	FORTRAN, PUBLIC, SUBROUTINE
*$ Declarations:
*	logical*4 function FPro_SincInfo(sinc, height, step, halfwidth,
*					shoulder)
*	integer*4	sinc
*	real*8		height, step, halfwidth, shoulder
*$ Inputs:
*	sinc		FORTRAN pointer to a sinc function object.
*$ Outputs:
*	height		peak y-value of function (at origin).
*	step		x-step between origin/zeroes of function.
*	halfwidth	x-value above which y-value is zero.
*	shoulder	x-width of truncation, modeled as a half-sinusoid.
*$ Returns:
*	.TRUE. if this is a sinc function; .FALSE. otherwise.
*$ Detailed_description:
*	This routine returns information about a sinc function object.
*$ External_references:
*	Profile toolkit
*$ Side_effects:
*	none
*$ Examples:
*	This snippet of code of code prints the step size in a sinc function.
*
*	integer*4	sinc
*	real*8		height, step, halfwidth, shoulder
*
*	call FPro_SincInfo(sinc, height, step, halfwidth, shoulder)
*	write(*,*) 'Step size = ', step
*$ Error_handling:
*	Profile toolkit error handling is in effect.
*
*	Conditions raised:
*	PRO_CLASS_ERROR		if this is not a software function object.
*	FORTRAN_POINTER_ERROR	if sinc is not a valid FORTRAN object pointer.
*$ Limitations:
*	none
*$ Author_and_institution:
*	Mark R. Showalter
*	NASA/Ames Research Center
*$ Version_and_date:
*	1.0: March 1998
*$ Change_history:
*	none
*******************************************************************************/

RL_INT4	FORTRAN_NAME(fpro_sincinfo) (sinc, height, step, halfwidth, shoulder)
RL_INT4	*sinc;
RL_FLT8	*height, *step, *halfwidth, *shoulder;
{
RL_VOID *ptr;
RL_BOOL	status;

    /* Look up sinc function pointer */
    ptr = FORT_GetPointer(*sinc);
    if (ptr == NULL) return 0;

    /* Return info */
    status = Pro_SincInfo(ptr, height, step, halfwidth, shoulder);
    return status ? FTRUE:FFALSE;
}

/*******************************************************************************