/* composit.c
********************************************************************************
* composit.c -- Routines for composite function objects
*
* User routines:
*	Pro_CompFunc(inner,outer)	creates a new function object that
*					returns a function of a function, i.e.
*					y = outer(inner(x)).
*
* Mark Showalter & Neil Heather, PDS Ring-Moon Systems Node, February 1998
*******************************************************************************/
#include <stdio.h>
#include <string.h>
#include "profile.h"
#include "fortran.h"

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

typedef struct ZPRO_COMPFUNC_STRUCT {
    XPRO_CLASS	class;
    PRO_OBJECT	*inner, *outer;
    RL_FLT8	outer_x1, outer_x2;
} ZPRO_COMPFUNC;

/********************
 * Static variables *
 ********************/

static XPRO_CLASS compfunc_class = {XPRO_COMPFUNC_CLASS, "composite", NULL};

/********************************
 * Internal function prototypes *
 ********************************/

static RL_FLT8 ZPro_EvalComp  RL_PROTO((RL_VOID *pointer, RL_FLT8 x));
static void    ZPro_FreeComp  RL_PROTO((RL_VOID *pointer));
static void    ZPro_PrintComp RL_PROTO((RL_VOID *pointer));

/*
********************************************************************************
* EXPORTED USER ROUTINES
********************************************************************************
*$ Component_name:
*	Pro_CompFunc (composit.c)
*$ Abstract:
*	This routine generates a new composite function object involving one
*	function of another.  The result returned is outer(inner(x)).
*$ Keywords:
*	PROFILE, FUNCTION
*	C, PUBLIC, SUBROUTINE
*$ Declarations:
*	PRO_OBJECT	*Pro_CompFunc(inner, outer)
*	PRO_OBJECT	*inner, *outer;
*$ Inputs:
*	inner		pointer to the inner function object (evaluated first).
*	outer		pointer to the outer function object (evaluated second).
*$ Outputs:
*	none
*$ Returns:
*	pointer to a new function object, or NULL on non-fatal error.
*$ Detailed_description:
*	This routine generates a new composite function object involving one
*	function of another.  The result returned is outer(inner(x)).
*$ External_references:
*	Profile toolkit
*$ Side_effects:
*	Memory is allocated.  The new object retains links both functions.
*$ Examples:
*	Suppose sqroot is a function object that returns the square root of x;
*	suppose square is a function object that returns the square of x.
*
*	func = Pro_CompFunc(square, sqroot);
*
*	Then Pro_FuncValue(func, 2.) returns 2.;
*	     Pro_FuncValue(func, 1.) returns 1.;
*	     Pro_FuncValue(func, 0.) returns 0.;
*	     Pro_FuncValue(func, -1.) returns 1.;
*$ Error_handling:
*	Profile library error handling is in effect.
*
*	Conditions raised:
*	PRO_CLASS_ERROR		if either argument is NULL or not a function.
*	PRO_COORD_MISMATCH	if the coordinate names do not match.
*	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_CompFunc(inner, outer)
PRO_OBJECT	*inner, *outer;
{
PRO_OBJECT	*new;
ZPRO_COMPFUNC	*compfunc;
RL_FLT8		x1, x2, size;

    /* Determine domain and validate inner function */
    size = Pro_ObjectDomain(inner, &x1,  &x2);
    if (size == 0.) return NULL;

    /* Test coordinate names */
    if (strcmp(Pro_ObjectName(outer,1), Pro_ObjectName(inner,2)) != 0) {
	XPro_CoordMismatch("composite function", inner, outer,
		Pro_ObjectName(inner,2), Pro_ObjectName(outer,1));
    }

    /* Create and initialize data structure */
    compfunc = (ZPRO_COMPFUNC *) XRL_Malloc(sizeof(ZPRO_COMPFUNC));
    if (compfunc == NULL) return NULL;

    compfunc->class = compfunc_class;
    compfunc->inner = inner;
    compfunc->outer = outer;
    (void) Pro_ObjectDomain(outer, &(compfunc->outer_x1),
                                   &(compfunc->outer_x2));

    /* Create object */
    new = XPro_MakeFunc(x1, x2, ZPro_EvalComp, ZPro_FreeComp, ZPro_PrintComp,
                        (RL_VOID *) compfunc);

    /* Enslave old functions and check for errors */
    if (new != NULL && (XPro_EnslaveObject(new,inner) ||
                        XPro_EnslaveObject(new,outer))) {
	Pro_FreeObject(new);
	new = NULL;
    }

    /* Transfer coordinate names to new object */
    Pro_RenameObject(new, 1, Pro_ObjectName(inner,1));
    Pro_RenameObject(new, 2, Pro_ObjectName(outer,2));

    return new;
}

/*
********************************************************************************
* INTERNAL ROUTINES
********************************************************************************
* ZPro_EvalComp(pointer, x)
*
* This internal routine calculates the value of a composite function.
*
* Inputs:
*	pointer		pointer to the ZPRO_COMPFUNC data structure.
*	x		location at which to evaluate composite function.
*
* Return:		value of composite function at x.
*
* Errors:
*	PRO_EVALUATION_FAILURE	if result of inner function falls outside
*				domain of outer function.
*******************************************************************************/

static RL_FLT8	ZPro_EvalComp(pointer, x)
RL_VOID		*pointer;
RL_FLT8		x;
{
ZPRO_COMPFUNC	*compfunc;
RL_FLT8		y1, y2;

    compfunc = (ZPRO_COMPFUNC *) pointer;

    /* Evaluate inner function */
    y1 = Pro_FuncValue(compfunc->inner, x);

    /* Check range for valid input to outer function */
    if (y1 < compfunc->outer_x1 || y1 > compfunc->outer_x2) {
	(void) sprintf(xpro_message, "function evaluation failure\n\
intermediate result falls outside domain in composite function\n\
f(%#g) = %#g in inner function \"%s\"\n\
domain = [%#g, %#g] in outer function \"%s\"",
		x, y1, Pro_ObjectName(compfunc->inner,0),
		compfunc->outer_x1, compfunc->outer_x2,
		Pro_ObjectName(compfunc->outer,0));

	RL_RaiseError("PRO_EVALUATION_FAILURE", xpro_message);

	if (y1 < compfunc->outer_x1) y1 = compfunc->outer_x1;
	else                         y1 = compfunc->outer_x2;
    }

    y2 = Pro_FuncValue(compfunc->outer, y1);

    return y2;
}

/*
********************************************************************************
* ZPro_FreeComp(pointer)
*
* This internal routine frees the memory used by a ZPRO_COMPFUNC structure.
*
* Inputs:
*	pointer		pointer to the ZPRO_COMPFUNC data structure.
*******************************************************************************/

static void	ZPro_FreeComp(pointer)
void		*pointer;
{
    XRL_Free(pointer);
}

/*
********************************************************************************
* ZPro_PrintComp(pointer)
*
* This internal routine prints out information about a composite function.  It
* is called by Pro_PrintFunc() and is used mainly for debugging.
*
* Input:
*	pointer		pointer to the ZPRO_COMPFUNC structure.
*******************************************************************************/

static void	ZPro_PrintComp(pointer)
RL_VOID        	*pointer;
{
ZPRO_COMPFUNC	*compfunc;

    compfunc = (ZPRO_COMPFUNC *) pointer;

    /* Make sure object is not NULL */
    if (compfunc == NULL) {
	printf("PRINT ERROR: Composite function pointer is NULL\n");
	return;
    }

    /* Make sure object is a composite function */
    if (compfunc->class.id != XPRO_COMPFUNC_CLASS) {
	printf("PRINT ERROR: Object is not a composite function\n");
	return;
    }

    /* Print object info... */
    printf("\nComposite function parameters...\n");
    printf("inner function = "); XPro_PrintInfo(compfunc->inner);
    printf("outer function = "); XPro_PrintInfo(compfunc->outer);

    return;
}

/*
********************************************************************************
* FORTRAN INTERFACE ROUTINES
********************************************************************************
*$ Component_name:
*	FPro_CompFunc (composit.c)
*$ Abstract:
*	This routine generates a new composite function object involving one
*	function of another.  The result returned is outer(inner(x)).
*$ Keywords:
*	PROFILE, FUNCTION
*	FORTRAN, PUBLIC, SUBROUTINE
*$ Declarations:
*	integer*4 function FPro_CompFunc(inner, outer)
*	integer*4	inner, outer
*$ Inputs:
*	inner		FORTRAN pointer to inner function (evaluated first).
*	outer		FORTRAN pointer to outer function (evaluated second).
*$ Outputs:
*	none
*$ Returns:
*	FORTRAN pointer to a new function object, or 0 on non-fatal error.
*$ Detailed_description:
*	This routine generates a new composite function object involving one
*	function of another.  The result returned is outer(inner(x)).
*$ External_references:
*	Profile toolkit
*$ Side_effects:
*	Memory is allocated.  The new object retains links both functions.
*$ Examples:
*	Suppose sqroot is a function object that returns the square root of x;
*	suppose square is a function object that returns the square of x.
*
*	func = FPro_CompFunc(square, sqroot)
*
*	Then FPro_FuncValue(func, 2.d0) returns 2.d0
*	     FPro_FuncValue(func, 1.d0) returns 1.d0
*	     FPro_FuncValue(func, 0.d0) returns 0.d0
*	     FPro_FuncValue(func, -1.d0) returns 1.d0
*$ Error_handling:
*	Profile library error handling is in effect.
*
*	Conditions raised:
*	PRO_CLASS_ERROR		if either argument is NULL or not a function.
*	PRO_COORD_MISMATCH	if the coordinate names do not match.
*	RL_MEMORY_ERROR		on memory allocation error.
*	FORTRAN_POINTER_ERROR	if object 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_compfunc) (inner, outer)
RL_INT4 *inner, *outer;
{
RL_VOID *ptr1, *ptr2, *ptr3;
RL_INT4 index;

    /* Look up function pointers */
    ptr1 = FORT_GetPointer(*inner);
    if (ptr1 == NULL) return 0;

    ptr2 = FORT_GetPointer(*outer);
    if (ptr2 == NULL) return 0;

    /* Call function */
    ptr3 = (RL_VOID *) Pro_CompFunc((PRO_OBJECT *) ptr1, (PRO_OBJECT *) ptr2);
    if (ptr3 == NULL) return 0;

    /* Save new pointer */
    index = FORT_AddPointer(ptr3);
    if (index == 0) Pro_FreeObject((PRO_OBJECT *) ptr3);

    return index;
}

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