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

  File name:  jemx_shared_imod.c
  Version:    see below
  Component:  
  Developer:  Peter Kretschmar, ISDC
              Peter.Kretschmar@obs.unige.ch
  
  Purpose:

  Utility functions to access JEM-X Instrument Model (IMOD) data.

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

/*
 * Revision 0.2.0  2004/01/23
 * First draft
 *
 */
#include "jemx_shared.h"




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

Function: JEMXdetectorGetMap

Description:
      
   Find the detector map for the given JEM-X instrument in the instrument
   model group and read in the data

Error codes:

0 == ISDC_OK   ==> function executed with no errors
< 0 (various)  ==> ISDC system function error code, passed on

Parameter arguments:
name           range         I/O     description

instModGrp     !=NULL       input   handle to instrument model group
jemxInstr      JMX1 or JMX2 input   instrument number 
detMap         (array)      output  2D array with DETE-MOD values
onAxisArea      >0.0        output  active area for an on-axis source
chatter        0 to 5       input   verbosity control variable
callingStatus  any          input   error status on input


*****************************************************************************/
int JEMXdetectorGetMap(dal_element     *instModGrp,
		       Instrument       jemxInstr,
		       unsigned short   detMap[J_NUM_X_PIXELS][J_NUM_Y_PIXELS],
		       double          *onAxisArea,
		       int              chatter,
		       int              callingStatus)
{
    /*
    *  Local variable definitions
    */
    int             status = ISDC_OK;       /* local status flag */
    dal_element    *detModPtr = NULL;       /* handle to JMXi-DETE-MOD */
    void           *buffer = NULL;
    long            dimAxes[2]={J_NUM_X_PIXELS,J_NUM_Y_PIXELS};
    long            ix=0, iy=0;
    char            detModName[14];

    /* exit function immediately if status is not OK on entry */
    if(callingStatus != ISDC_OK) return(callingStatus);
    
    switch(jemxInstr) {
    case JMX1:
	strcpy(detModName,"JMX1");
	break;
    case JMX2:
	strcpy(detModName,"JMX1");
	break;
    default:
	status = JEMX_SHARED_ILLEGAL_DETECTOR_NUMBER; 
	RILlogMessage(NULL,Error_1,
		      "JEM-X instrument number %i impossible!",jemxInstr);
	break;
    }
    strcat(detModName,J_DETECTOR_DS);
    
    status = JEMXarrayGet(instModGrp,detModName,DAL_USHORT,
			  2,dimAxes,&detModPtr,&buffer,
			  chatter,status);
    if (ISDC_OK!=status) {
	return(status);
    }

    /* Map 1D transfer buffer to 2D array */
    for (ix=0; ix<J_NUM_X_PIXELS; ix++) {
	for (iy=0; iy<J_NUM_Y_PIXELS; iy++) {
	    detMap[ix][iy] = 
		*( (unsigned short *)buffer + ix + iy * J_NUM_X_PIXELS );
	}
    }

    /*
     * Read area attribute, 
     * if unknown (older IC data) set to default constant from jemx.h
     */
    status = DALattributeGetReal(detModPtr,"ONAXAREA",onAxisArea,
				 NULL,NULL,status);
    if (CF_KEY_NO_EXIST==status) {
	status = ISDC_OK;
	*onAxisArea = J_IDEAL_ILLUM;
	if (J_CHATTY_NORMAL<=chatter) {
	    RILlogMessage(NULL,Warning_2,
			  "On axis area keyword not found in %s! "
			  "Will use default value %f",
			  detModName,*onAxisArea);
	}
    } else if (ISDC_OK!=status) {
	*onAxisArea = 0.0;
	RILlogMessage(NULL,Error_2,
		      "Error %d: reading On axis area keyword from %s!",
		      status,detModName);
	return(status);
    } else if (J_CHATTY_VERY_VERBOSE<=chatter) {
	RILlogMessage(NULL,Log_0,
		      "On axis area %lf read from %s",*onAxisArea,detModName);
    }

    /* Free used memory */
    if (NULL!=buffer) { DALfreeDataBuffer(buffer,ISDC_OK); }

    return(status);
}



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

Function: JEMXdetectorGetCorXY

Description:
      
   Read in the two tables with the corrected pixel positions along
   the X and Y axis of the detector.

Error codes:

0 == ISDC_OK   ==> function executed with no errors
< 0 (various)  ==> ISDC system function error code, passed on

Parameter arguments:
name           range         I/O     description

instModGrp     !=NULL       input   handle to instrument model group
jemxInstr      JMX1 or JMX2 input   instrument number 
detCorX        (array)      output  2D array with corrected X positions
detCorY        (array)      output  2D array with corrected Y positions
chatter        0 to 5       input   verbosity control variable
callingStatus  any          input   error status on input


*****************************************************************************/
int JEMXdetectorGetCorXY(dal_element *instModGrp,
			 Instrument   jemxInstr,
			 double       detCorX[J_NUM_X_PIXELS][J_NUM_Y_PIXELS],
			 double       detCorY[J_NUM_X_PIXELS][J_NUM_Y_PIXELS],
			 int          chatter,
			 int          callingStatus)
{
    /*
    *  Local variable definitions
    */
    int             status = ISDC_OK;       /* local status flag */
    dal_element    *detCorPtr = NULL;       /* handle to JMXi-DETE-MOD */
    void           *buffer = NULL;
    long            dimAxes[2]={J_NUM_X_PIXELS,J_NUM_Y_PIXELS};
    long            ix=0, iy=0;
    char            detCorName[14];

    /* exit function immediately if status is not OK on entry */
    if(callingStatus != ISDC_OK) return(callingStatus);
    
    switch(jemxInstr) {
    case JMX1:
	strcpy(detCorName,"JMX1");
	break;
    case JMX2:
	strcpy(detCorName,"JMX1");
	break;
    default:
	status = JEMX_SHARED_ILLEGAL_DETECTOR_NUMBER; 
	RILlogMessage(NULL,Error_1,
		      "JEM-X instrument number %i impossible!",jemxInstr);
	break;
    }

    /* Get CORX table */
    strcpy(detCorName+4,J_CORX_DS);
    
    status = JEMXarrayGet(instModGrp,detCorName,DAL_DOUBLE,
			  2,dimAxes,&detCorPtr,&buffer,
			  chatter,status);
    if (ISDC_OK!=status) {
	return(status);
    }
   
    /* Map 1D transfer buffer to 2D array */
    for (ix=0; ix<J_NUM_X_PIXELS; ix++) {
	for (iy=0; iy<J_NUM_Y_PIXELS; iy++) {
	    detCorX[ix][iy] = 
		*( (double *)buffer + ix + iy * J_NUM_X_PIXELS );
	}
    }

    /* Get CORY table */
    strcpy(detCorName+4,J_CORY_DS);
    detCorPtr = NULL;

    status = JEMXarrayGet(instModGrp,detCorName,DAL_DOUBLE,
			  2,dimAxes,&detCorPtr,&buffer,
			  chatter,status);
    if (ISDC_OK!=status) {
	return(status);
    }
   
    /* Map 1D transfer buffer to 2D array */
    for (ix=0; ix<J_NUM_X_PIXELS; ix++) {
	for (iy=0; iy<J_NUM_Y_PIXELS; iy++) {
	    detCorY[ix][iy] = 
		*( (double *)buffer + ix + iy * J_NUM_X_PIXELS );
	}
    }

    /* Free used memory */
    if (NULL!=buffer) { DALfreeDataBuffer(buffer,ISDC_OK); }

    return(status);
}

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

Function: JEMXdetectorPixelAreas

Description:
      
   Calculate the lower left corners and the areas of the (corrected)
   pixels.
   

Error codes:

0 == ISDC_OK   ==> function executed with no errors
< 0 (various)  ==> ISDC system function error code, passed on

Parameter arguments:
name           range         I/O     description

detCorX        (array)      input   2D array with corrected X positions
detCorY        (array)      input   2D array with corrected Y positions
pixelArea      (array)      output  2D array with effective sizes of pixels
pixelLowLeftX  (array)      output  2D array with pixel lower left X pos.
pixelLowLeftY  (array)      output  2D array with pixel lower left Y pos.
chatter        0 to 5       input   verbosity control variable
callingStatus  any          input   error status on input


*****************************************************************************/
int JEMXdetectorPixelAreas(double  detCorX[J_NUM_X_PIXELS][J_NUM_Y_PIXELS],
			   double  detCorY[J_NUM_X_PIXELS][J_NUM_Y_PIXELS],
			   double  pixelArea[J_NUM_X_PIXELS][J_NUM_Y_PIXELS],
			   double  pixelLowLeftX[J_NUM_X_PIXELS][J_NUM_Y_PIXELS],
			   double  pixelLowLeftY[J_NUM_X_PIXELS][J_NUM_Y_PIXELS],
			   int     chatter,
			   int     callingStatus)
{
    /*
    *  Local variable definitions
    */
    int             status = ISDC_OK;       /* local status flag */
    /*
     *
     *  Pixel (x,y) with corners:
     *
     *
     *  (x-1,y+1)         (x,y+1)            (x+1,y+1)
     *
     *        (ulCorner)     (urCorner)
     *
     *  (x-1,y)            (x,y)             (x+1,y)
     *
     *        (llCorner)    (lrCorner)
     *
     *  (x-1,y-1)         (x,y-1)            (x+1,y-1)
     *
     */
    struct pixelCorner { 
	double x; 
	double y; 
    } llCorner, ulCorner, lrCorner, urCorner;

    double height    = 0.0;
    double width     = 0.0;
    double rawRadius = 0.0;
    double sumArea   = 0.0;

    long ix=0, iy=0;

    /* exit function immediately if status is not OK on entry */
    if(callingStatus != ISDC_OK) return(callingStatus);
    
    /* 
     *  Loop over all detector pixels. 
     *  Set corners as mean of surrounding pixel coordinates.
     *  Calculate pixel area assuming a rectangle
     */
    for ( ix=1 ; ix<J_NUM_X_PIXELS-1  ; ix++){
	for ( iy=1 ; iy<J_NUM_X_PIXELS-1  ; iy++){
	    llCorner.x = 0.25*(detCorX[ix-1][iy-1] +
			       detCorX[ix-1][iy] +
			       detCorX[ix][iy] +
			       detCorX[ix][iy-1]);
	    llCorner.y = 0.25*(detCorY[ix-1][iy-1] +
			       detCorY[ix-1][iy] +
			       detCorY[ix][iy] +
			       detCorY[ix][iy-1]);
	    ulCorner.x = 0.25*(detCorX[ix-1][iy+1] +
			       detCorX[ix][iy+1] +
			       detCorX[ix][iy] +
			       detCorX[ix-1][iy]);
	    ulCorner.y = 0.25*(detCorY[ix-1][iy+1] +
			       detCorY[ix][iy+1] +
			       detCorY[ix][iy] +
			       detCorY[ix-1][iy]);
	    urCorner.x = 0.25*(detCorX[ix+1][iy+1] +
			       detCorX[ix][iy+1] +
			       detCorX[ix][iy] +
			       detCorX[ix+1][iy]);
	    urCorner.y = 0.25*(detCorY[ix+1][iy+1] +
			       detCorY[ix][iy+1] +
			       detCorY[ix][iy] +
			       detCorY[ix+1][iy]);
	    lrCorner.x = 0.25*(detCorX[ix+1][iy] +
			       detCorX[ix+1][iy-1] +
			       detCorX[ix][iy] +
			       detCorX[ix][iy-1]);
	    lrCorner.y = 0.25*(detCorY[ix+1][iy] +
			       detCorY[ix+1][iy-1] +
			       detCorY[ix][iy] +
			       detCorY[ix][iy-1]);

	    pixelLowLeftX[ix][iy] = llCorner.x;
	    pixelLowLeftY[ix][iy] = llCorner.y;


	    /* 
	     * Estimate pixel area by rectangle approximation 
	     * converting from the mm coordinates to cm.
	     */
	    
	    height = 0.1 * (0.5*(ulCorner.y+urCorner.y) - 
			    0.5*(llCorner.y+lrCorner.y));
	    width  = 0.1 * (0.5*(lrCorner.x+urCorner.x) - 
			    0.5*(llCorner.x+ulCorner.x));

	    /* Get pixel position off center */
	    rawRadius = sqrt(
		(double) ((ix-JEMX_RAW_CENTER_X)*(ix-JEMX_RAW_CENTER_X)) +
		(double) ((iy-JEMX_RAW_CENTER_Y)*(iy-JEMX_RAW_CENTER_Y)) );

	    if(rawRadius < JEMX_MAX_USEFUL_RADIUS){
		pixelArea[ix][iy] = height*width;
		sumArea += pixelArea[ix][iy];
	    }else{
		pixelArea[ix][iy] = 0.01; /* 1mm x 1mm */
	    }

	    
	} /* end loop iy */
    } /* end loop ix */

    if (J_CHATTY_VERBOSE<chatter) {
	RILlogMessage(NULL,Log_2,
		      "The total area within a raw radius of %f is %f.",
		      JEMX_MAX_USEFUL_RADIUS,sumArea);
    }

    return(status);
}


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

Function: JEMXdetectorGetData

Description:
      
   Read in the relevant information for the JEM-X detector.

Error codes:

0 == ISDC_OK   ==> function executed with no errors
< 0 (various)  ==> ISDC system function error code, passed on

Parameter arguments:
name           range         I/O     description

instModGrp     !=NULL       input   handle to instrument model group
jemxInstr      JMX1 or JMX2 input   instrument number 
detector       (struct)     output  important detector data
chatter        0 to 5       input   verbosity control variable
callingStatus  any          input   error status on input


*****************************************************************************/
int JEMXdetectorGetData(dal_element      *instModGrp,
			Instrument        jemxInstr,
			JEMXdetectorData *detector,
			int               chatter,
			int               callingStatus)
{
    /*
    *  Local variable definitions
    */
    int              status = ISDC_OK;       /* local status flag */

    /* exit function immediately if status is not OK on entry */
    if(callingStatus != ISDC_OK) return(callingStatus);

    /* read in pixel status information */
    status =  JEMXdetectorGetMap(instModGrp,jemxInstr,
				 detector->detMap,&(detector->area),
				 chatter,status);
    if (ISDC_OK!=status) {
	RILlogMessage(NULL,Error_2,
		      "Error %d reading detector map",status);
	return(status);
    }
    
    /* read in the corrected pixel positions */
    status =  JEMXdetectorGetCorXY(instModGrp,jemxInstr,
				   detector->CorX,detector->CorY,
				   chatter,status);
    if (ISDC_OK!=status) {
	RILlogMessage(NULL,Error_2,
		      "Error %d reading CORX/CORY tables",status);
	return(status);
    }
    
    /* derive the areas per pixel */
    status = JEMXdetectorPixelAreas(detector->CorX,detector->CorY,
				    detector->pixelArea,
				    detector->pixelLowerLeftX,
				    detector->pixelLowerLeftY,
				    chatter,status);

    return(status);
} 


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

Function: JEMXdetectorCalcUsedArea

Description:
      
   Calculate the overall used area within the detector based on the given
   radius limit and excluding bad pixels.

Error codes:

0 == ISDC_OK   ==> function executed with no errors
< 0 (various)  ==> ISDC system function error code, passed on

Parameter arguments:
name           range         I/O      description

detector       (struct)     input     detector description
radiusLimit    >=0.0        input     radius in mm up to which events are used
excludeFlags   >=0          input     holds bit flags excluding pixels
area           >=0.0        output    effectively used area in cm^2
usedMap        (array)      output    
chatter        0 to 5       input     verbosity control variable
callingStatus  any          input     error status on input


*****************************************************************************/
int JEMXdetectorCalcOnAxisArea(JEMXdetectorData  detector,
			       double            radiusLimit,
			       unsigned short    excludeFlags,
			       double           *area,
			       unsigned short    
			       usedMap[J_NUM_X_PIXELS][J_NUM_Y_PIXELS],
			       int               chatter,
			       int               callingStatus)
{
    /*
    *  Local variable definitions
    */
    int             status = ISDC_OK;       /* local status flag */
    unsigned short  excludePixel = 0;
    long            ix=0, iy=0, usedPixels=0;
    double          radius=0.0;

    /* exit function immediately if status is not OK on entry */
    if(callingStatus != ISDC_OK) return(callingStatus);

     /* 
     *  Loop over all detector pixels. 
     *  Exclude those with certain bits set.
     */
    *area = 0.0;
    usedPixels = 0;

    for ( ix=1 ; ix<J_NUM_X_PIXELS-1  ; ix++){
	for ( iy=1 ; iy<J_NUM_X_PIXELS-1  ; iy++){
	
	    usedMap[ix][iy] = 0;

	    /* Radius from center of detector */
	    radius = sqrt( (detector.CorX[ix][iy]*detector.CorX[ix][iy]) +
			   (detector.CorY[ix][iy]*detector.CorY[ix][iy]) );

	    if (radius <= radiusLimit) {
		/* Compare bits of excludeFlages with detector map data */
		excludePixel = ( excludeFlags & detector.detMap[ix][iy] );
		if (excludePixel) {
		    if (J_CHATTY_MAXIMAL <= chatter) {
			RILlogMessage(NULL,Log_0,
				      "Excluding pixel (%ld,%ld) with value %d",
				      ix+1,iy+1,detector.detMap[ix][iy]);
		    }
		} else {
		    *area += detector.pixelArea[ix][iy];
		    usedMap[ix][iy] = detector.detMap[ix][iy];
		    usedPixels++;
		}
	    }

	} /* end loop iy */
    } /* end loop ix */
    if (J_CHATTY_VERBOSE <= chatter) {
	RILlogMessage(NULL,Log_1,
		      "Using %ld detector pixels for a total detector area of %f",
		      usedPixels,*area);
    }

    return(status);
}