//#############################################################################
// ASL IPS Compute Ice Draft Routines
//
// (c) CopyRight ASL Environmental Sciences Inc. 1997-2001
//
//	File:			calcAvgProp.c
//
//	Purpose:		Calculate the average sound speed and density using a CTD profile
//				for the calculation of ice drafts.
//
//	Development:		National Instruments Labwindows/CVI 5.0 and higher
//					
//	Date:			March 20, 2001
//
// 	Programmer:		Dave Billenness
//
//  	Modified 		April 11, 2001 - Changed definitions of date to short ints
//				Changed to only one include file "IceDraft.h" "IceDraft_def.h
//				removed.
//				April 27, 2001 - Function to calc avg SOS and density separated
//				from main function IceDraft.c
//******************************************************************************
//#include <ansi_c.h>
#include <stdio.h>
#include <math.h>
#include "IceDraft.h"


//#############################################################################
// Section:  Ice Draft Programs
//
// Function: int calcAvgProp(double Patm, int N, double Conductivity[], 
//			double Temperature[], double Pressure[], double *AverageSoundSpeed, double *AverageDensity)
//
// Purpose: Compute average sound speed and average density given:
//			Patm - estimated value of the atmospheric pressure (in bars) at the time of the IPS measurement
//			N - number of elements in CTD array
//			C - pointer to first element of conductivity array (C(S,T,0)/C(35,T,0))
//			T - pointer to first element of temperature array (degress Celsius)
//			P - pointer to first element of pressure array (absolute pressure in bars)
//			AverageSoundSpeed - computed value of the average sound speed (in m/s)
//			AverageDensity - computed value of the average density (in kg/m^3)
//
// returns: an int ErrorCode as defined in IceDraft.h
//
//#############################################################################
int calcAvgProp(double Patm, int N, double Conductivity[], double Temperature[], double Pressure[], double *AverageSoundSpeed, double *AverageDensity)
{
		int 			i = 0, Error = 0;
		static double	Salinity[ICEDRAFT_MAXARRAYSIZE], ScaledPressure[ICEDRAFT_MAXARRAYSIZE], SoundSpeed[ICEDRAFT_MAXARRAYSIZE], Density[ICEDRAFT_MAXARRAYSIZE];
		double			MinPressure = 0., MaxPressure = 0., AvgSoundSpeed = 0., AvgDensity = 0.;
		
		
		// Now start to look at CTD profile to compute an average sound speed and density.  
		// Convert CTD Pressure in bars to decibars using CONVERT_PRESSURE
		// Also convert pressure to gauge by removing Patm
		for(i=0; i<N; i++)
		{
			ScaledPressure[i] = (Pressure[i] - Patm)*CONVERT_PRESSURE;
		}	
		
		// Compute Salinity from conductivity ratio, pressure and temperature using pss78 code
		// then compute Sound Speed and Density. 
		// At the same time find the min and max pressure for setting up the pressure bins.
		MaxPressure = ScaledPressure[0];
		MinPressure = ScaledPressure[0];
		for(i=0; i<N; i++)
		{
			if(ScaledPressure[i] > MaxPressure)
				MaxPressure = ScaledPressure[i];
			if(ScaledPressure[i] < MinPressure)
				MinPressure = ScaledPressure[i];
				
			// Compute Salinity	
			Error = pss78(&Salinity[i], ScaledPressure[i], Temperature[i], Conductivity[i]);
			
			// Check for a valid salinity - if salinity is negative then function State will crash
			if(Error < 0) 
				return Error;
			if(Salinity[i] > 50)
				return ERROR_CTD_INVALID_SALINITY;
			
			// Compute SoundSpeed				
			Error = sound(&SoundSpeed[i], ScaledPressure[i], Temperature[i], Salinity[i]);
			
			// Check for a valid SoundSpeed
			if(Error < 0)
				return Error;
		
			// compute Density	
			Error = state(&Density[i],1, ScaledPressure[i], Temperature[i], Salinity[i]);
			
			// Check for a valid Density
			if(Error < 0)
				return Error;
		} //for(i=0; i<N; i++)
		
		// Now go through the Pressure, Salinity, Sound Speed and Density arrays and
		// bin them based on the Pressure - calculate the min and max Pressure and
		// the number of depth bins and place Sal, SS and Density into bins
		Error = bins(MinPressure, MaxPressure, ScaledPressure, SoundSpeed, Density, N, &AvgSoundSpeed, &AvgDensity);
		if(Error != 0)
			return Error;
				
		
		// Check AvgDensity
		if(AvgDensity == 0)
			return ERROR_DENSITY_EQUAL_0;
		
		*AverageDensity = AvgDensity;
		*AverageSoundSpeed = AvgSoundSpeed;
		
		return Error;
}		



//#############################################################################
// Section:  IceDraft  Program
//
// Function: int bins(double min, double max, double Data[], double A[], double B[], int N, double *AvgA, double *AvgB)
//
// Purpose:  put soundspeed and density data into pressure bins then average, interpolating values with
//			no data in the bins
//#############################################################################
int bins(double min, double max, double Data[], double A[], double B[], int N, double *AvgA, double *AvgB)
{
	int				Error = 0;						// Error Code
	long int 		histonum[ICEDRAFT_MAXARRAYSIZE];// Array with number of values in each bin
	long int 		bin_num;						// Bin number for histogram
	long int	 	i = 0, j = 0, NumBins = 0;		// Loop counter, Number of Bins
	double 			temp;
	static double	Data1[ICEDRAFT_MAXARRAYSIZE], Data2[ICEDRAFT_MAXARRAYSIZE], Data3[ICEDRAFT_MAXARRAYSIZE];

	
	// Calculate the number of bins based on the bin width (global variable)
	NumBins = (ceil(max) - floor(min))/BINWIDTH + 1;

	// NumBins must be <= ICEDRAFT_MAXARRAYSIZE
	if(NumBins > ICEDRAFT_MAXARRAYSIZE)
		return ERROR_EXCEEDED_MAX_CTD_ARRAYSIZE;
	
	// Initialize the arrays to zero
	for (i=0; i<NumBins; i++)
	{
  		histonum[i] = 0;
  		Data1[i] = 0.;
  		Data2[i] = 0.;
  		Data3[i] = 0.;
  	}
  
	if(NumBins > 0)
	{
		// Bin the sample
		for (i=0; i<N; i++)
		{

    		// Map data to its appropriate bin
	    	for (bin_num=0; bin_num<=NumBins; bin_num++)
	    	{
      			if (Data[i] < ((double) bin_num+floor(min)))
      			{
					// Increment selected bin_number
	    			histonum[bin_num-1]++;
    				Data1[bin_num-1] = Data1[bin_num-1] + Data[i];
    				Data2[bin_num-1] = Data2[bin_num-1] + A[i];
    				Data3[bin_num-1] = Data3[bin_num-1] + B[i];
      				
      				break;
      			}
      		} //for (bin_num=0; bin_num<=NumBins; bin_num++)	
      		
		} //for (i=0; i<N; i++)

		// Average the values in the bins using the number of values
		// If a bin is zero then use the previous value
		for (bin_num=0; bin_num<NumBins; bin_num++)
		{
			if(Data1[bin_num] != 0)
			{
				Data1[bin_num] = Data1[bin_num]/histonum[bin_num];
				Data2[bin_num] = Data2[bin_num]/histonum[bin_num];
				Data3[bin_num] = Data3[bin_num]/histonum[bin_num];
				j++;
				
				// Now sum Data2 and Data3
				*AvgA = *AvgA + Data2[bin_num];
				*AvgB = *AvgB + Data3[bin_num];
			} //if(Data1[bin_num] != 0) 
		} //for (bin_num=0; bin_num<NumBins; bin_num++)
		
		// Now average Data2 and Data3 - this average is the average of the non-zero values
		*AvgA = *AvgA/j;
		*AvgB = *AvgB/j;
		
		// Interpolate 0 values (areas with no bin data)
		// and then come up with an average over the entire range.  If a 0 appears at
		// the beginning of end of the array, replace with AvgA or AvgB
		RemoveZeroes(Data2, *AvgA, 0, NumBins);
		RemoveZeroes(Data3, *AvgB, 0, NumBins);
		
		// Now get an average of all of the data
		*AvgA = 0;
		*AvgB = 0;
		for (bin_num=0; bin_num<NumBins; bin_num++)
		{
			*AvgA = *AvgA + Data2[bin_num];
			*AvgB = *AvgB + Data3[bin_num];
	
		} //for (bin_num=0; bin_num<NumBins; bin_num++)
		
		// Now average Data2 and Data3 - this average is the average of all values, including interpolated values
		*AvgA = *AvgA/NumBins;
		*AvgB = *AvgB/NumBins;
		
	}else
	{
		*AvgA = A[0];
		*AvgB = B[0];
	}
	

	return Error;
}



//#############################################################################
// Section:  IceDraft Program
//
// Function: void RemoveZeroes(float *DataPtr, float Mean, int Value2Interp, long int N)
//
// Purpose:  interpolate values=Value2Interp in an array.  Replace points at the
//			beginning and end with the mean
//#############################################################################
void RemoveZeroes(double *DataPtr, double Mean, int Value2Interp, long int N)
{
		
		int 	i, j, k, Pts2Interp;
		double 	FirstLocX, EndLocX, *FirstLocY, *EndLocY;
		double 	x[2], y[2];
		double	Interp_Val;
		
		// loop through the Data and remove zeroes
		for(i = 0; i < N; i++)
		{
			Pts2Interp = 0;
			
			if(*DataPtr == Value2Interp)
			{
			
				// This could be the first point.  If it is then we have no starting point
				// to interpolate to and the best we can do is replace it with the mean
				if(i == 0)
				{
					*DataPtr = Mean;
					DataPtr++;
					continue;
				
				} //if(i == 0)
			
				// remember this start good location
				FirstLocY = (DataPtr-1);
				FirstLocX = i-1;
				
				// we have found a zero - now loop until we find a non zero
				while(*DataPtr == Value2Interp)
				{
					// this might be the last point in the data - if the last set of 
					// points are zero
					if(i == N-1)
					{
						// set the last point to the previous value of the burst
						*DataPtr = *FirstLocY; //Mean;
						break;
					}	
					DataPtr++;
					Pts2Interp++;
					i++;
				} //while(*DataPtr == Value2Interp)
				
				// remember this end bad location
				EndLocY = DataPtr;
				EndLocX = i;
				
				// we know the location of the start and end points now we can interpolate
				for(j = 0; j < Pts2Interp; j++)
				{
					x[0] = FirstLocX;
					x[1] = EndLocX;
					y[0] = *FirstLocY;
					y[1] = *EndLocY;
					Interp_Val = ((FirstLocX+j+1-x[0])/(x[1] - x[0])) * (y[1] - y[0]) + y[0];
					
					// place the value in the Data
					*(FirstLocY+j+1) = Interp_Val;
				}  //for(j = 0; j < Pts2Interp; j++)
				
			} //if(*DataPtr == Value2Interp)	
			
			DataPtr++;
		
		}  //for(i=0; i<N; i++)
		
		return;
}