/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*                                                                        */
/*  Incluse file:     ArrayCalc.c                                         */
/*                                                                        */
/*  Main file:        BESmethod.c                                         */
/*                                                                        */
/*  Author:   Dr.A.B.Giles    RXTE PCA group, USRA, NASA / GSFC           */
/*                                                                        */
/*  Date:     June 21 1997                                                */
/*                                                                        */
/*  Version:  1.0                                                         */
/*                                                                        */
/*  Contact:  Mail:   Code 662, NASA / GSFC, Greenbelt, MD 20771, USA     */
/*            Tel:    301 286 6628                                        */
/*            FAX:    301 286 1684                                        */
/*            E-mail  barry@pcasun3.gsfc.nasa.gov                         */
/*                                                                        */
/*                                                                        */
/*  Suite of routines to manipulate the various analysis arrays           */
/*                                                                        */
/*  void ProcessAnalysis(int BinHeight)                                   */
/*  void EvaluateBurst(int StatRow, int BurstHeight, long WindowSum)      */
/*  void ListAnalysisSummary(int KeyPressRequired, int ClrScrRequired)    */
/*  void ListAnalysis(int KeyPressRequired, int ClrScrRequired)           */
/*  void ListResults(int KeyPressRequired, int ClrScrRequired)            */
/*  void ListNormalisation(int KeyPressRequired, int ClrScrRequired)      */
/*  void ListScale(int KeyPressRequired, int ClrScrRequired)              */
/*  void ListAll()                                                        */
/*  void ClearResults()                                                   */
/*  void ClearAnalysis()                                                  */
/*  void FileWriteResults(int num)                                        */
/*  void SaveLightCurve(int num)                                          */
/*                                                                        */
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/



void ProcessAnalysis(BinHeight)
int BinHeight;
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Process new bin height value into the analysis data arrays  */
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	int j,jj,k1,k2,temp;
	int BurstPosn,BurstHeight;
	long WindowSum;
	float olddata,newdata;

	temp = AnalysisBin[0];
	olddata = Analysis[0][temp];
	Analysis[0][temp] = BinHeight;
	AnalysisMark[0] = AnalysisMark[0] + 1;
	AnalysisSum[0] = AnalysisSum[0] + BinHeight - olddata;
	AnalysisBin[0] = (temp + 1) & AnalysisWidth_1;
	AnalysisCen[0] = (AnalysisCen[0] + 1) & AnalysisWidth_1;
	if (AnalysisLive == 1) {
		BurstPosn = AnalysisCen[0];
		BurstHeight = Analysis[0][BurstPosn];
		WindowSum = AnalysisSum[0];
		jj = 0;
		EvaluateBurst(jj,BurstHeight,WindowSum);
	}
	for (j=0; j<ResultHeight_1; j++) {
		if (AnalysisMark[j] == 2) {
			AnalysisMark[j] = 0;
			temp = AnalysisBin[j];
			k1 = (temp - 2) & AnalysisWidth_1;
			k2 = (temp - 1) & AnalysisWidth_1;
			newdata = Analysis[j][k1] + Analysis[j][k2];
			jj = j + 1;
			temp = AnalysisBin[jj];
			olddata = Analysis[jj][temp];
			Analysis[jj][temp] = newdata;
			AnalysisMark[jj] = AnalysisMark[jj] + 1;
			AnalysisSum[jj] = AnalysisSum[jj] + newdata - olddata;
			AnalysisBin[jj] = (temp + 1) & AnalysisWidth_1;
			AnalysisCen[jj] = (AnalysisCen[jj]+1) & AnalysisWidth_1;
			if (AnalysisLive == 1) {
				BurstPosn = AnalysisCen[jj];
				BurstHeight = Analysis[jj][BurstPosn];
				WindowSum = AnalysisSum[jj];
				EvaluateBurst(jj,BurstHeight,WindowSum);
			}
		}
	}
	if (AnalysisMark[ResultHeight_1] == AnalysisWidth) {
		AnalysisLive = 1;
		AnalysisMark[ResultHeight_1] = 0;
	}
}



void EvaluateBurst(StatRow, BurstHeight, WindowSum)
int StatRow, BurstHeight;
long WindowSum;
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Evaluate burst height (bin value) against window sum (mean)  */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	int i,j,level,stat,lev;
	double WindowAverage;
	float TableMean;
	char ch;
	FILE *zp;

/*  Increment try array  */
	AnalysisTry[StatRow] = AnalysisTry[StatRow] + 1;
/*  Compute average over window  */
	WindowAverage = (WindowSum * 1.0) / AnalysisWidth;
/*  Constrain burst height to valid computed table probabilities  */
	if (BurstHeight >= MaxXrayHt) BurstHeight = MaxXrayHt - 1;

       if (BurstHeight >= StartOfTable) {
/*  Perform test - only increment the most significant probability bin  */
	level = -1;
	for (i=0; i<ResultWidth; i++) {
		TableMean = Comp[BurstHeight][i];
		if ((WindowAverage < TableMean) && (level < i)) level = i;
	}

/*  Level > -1 if at least one of the probability tests was positive  */
	if (level > -1) {

/*  Possible test output of all significant bursts  */
		Result[StatRow][level] = Result[StatRow][level] + 1;
		if (Test_Detections == 'Y') {
		  printf(" %6d  %6d  %6d \n",StatRow,BurstHeight,WindowSum);
		  for (i=0; i<ResultWidth; i++) {
			printf(" %7.4f",Comp[BurstHeight][i]);
		  }
		  printf("     %7.4f  %7d \n",WindowAverage,level);
		  ch = getch();
		}

/*  Individual significant burst details if above 'Example limit' setting  */
		if (Elog == 1) {
			stat = StatRow + 1;
			lev = level + 1;
			if (lev >= ExampleLimit) {
				NoExamples = NoExamples + 1;
				fprintf(eecon," No. =%5d",NoExamples);
				fprintf(eecon,"  Bin No. =%12.1f",SuperTime);
				fprintf(eecon,"  Row = %2d",stat);
				fprintf(eecon,"  Ht = %4d",BurstHeight);
				fprintf(eecon,"  Av = %7.4f",WindowAverage);
				fprintf(eecon,"  Col = %2d\n",lev);
			}
		}
	}
      }

/*  Possible test output for every visit to this routine  */
	if (Test_Evaluate == 'Y') {
		printf(" %6d  %6d  %6d \n",StatRow,BurstHeight,WindowSum);
		for (i=0; i<ResultWidth; i++) {
			printf(" %7.4f",Comp[BurstHeight][i]);
		}
		printf(" %7.4f  %7d \n",WindowAverage,level);
	}
}



void ListAnalysisSummary(KeyPressRequired, ClrScrRequired)
int KeyPressRequired, ClrScrRequired;
/*+++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  List the contents of the analysis summary array  */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	int i,j,k,tempWidth,BurstSize;
	float ff;
	char ch;

	if (ClrScrRequired == 1) ClearScreen();
	printf("Working Array Summary:-\n\n");
	printf("       Bin     Burst     Burst       Row       No.\n");
	printf("   Pointer   Pointer     Value       Sum      Trys\n");
	for (i=0; i<ResultHeight; i++) {
		printf("%10d",AnalysisBin[i]);
		printf("%10d",AnalysisCen[i]);
		j = AnalysisCen[i];
		BurstSize = Analysis[i][j];
		printf("%10d",BurstSize);
		ff = AnalysisSum[i];
		printf("%10.0f",ff);
		printf("%10d",AnalysisTry[i]);
		printf("\n");
	}
	printf("\r\n\n");
	if (KeyPressRequired == 1) ch = getch();
}



void ListAnalysis(KeyPressRequired, ClrScrRequired)
int KeyPressRequired, ClrScrRequired;
/*+++++++++++++++++++++++++++++++++++++++++++*/
/*  List the contents of the analysis array  */
/*+++++++++++++++++++++++++++++++++++++++++++*/
{
	int i,j,k,tempWidth;
	char ch;

	if (ClrScrRequired == 1) ClearScreen();
	if (AnalysisLive == 0) printf("Working Array:- (LOADING)\n\n");
	if (AnalysisLive == 1) printf("Working Array:- (LIVE)\n\n");
	tempWidth = AnalysisWidth;
	if (tempWidth > 12) tempWidth = 12;
	for (i=0; i<ResultHeight; i++) {
		for (j=0; j<tempWidth; j++) {
			k = Analysis[i][j];
			printf(" %5d",k);
		}
		printf("\n");
	}
	printf("\r\n\n");
	if (KeyPressRequired == 1) ch = getch();
}



void ListResults(KeyPressRequired, ClrScrRequired)
int KeyPressRequired, ClrScrRequired;
/*++++++++++++++++++++++++++++++++++++++++++*/
/*  List the contents of the results array  */
/*++++++++++++++++++++++++++++++++++++++++++*/
{
	int i,j,k;
	char ch;

	if (ClrScrRequired == 1) ClearScreen();
	if (AnalysisLive == 0) 
          printf("Result Array - No. of Identified Bursts:- (LOADING)\n\n");
	if (AnalysisLive == 1)
	  printf("Result Array - No. of Identified Bursts:- (LIVE)\n\n");
	for (i=0; i<ResultWidth; i++) printf(" %7.4f",Prob_Limit_Array[i]);
	printf("\n");
	for (i=0; i<ResultHeight; i++) {
		for (j=0; j<ResultWidth; j++) {
			k = Result[i][j];
			printf(" %7d",k);
		}
		printf("\n");
	}
	printf("\r\n\n");
	if (KeyPressRequired == 1) ch = getch();
}



void ListNormalisation(KeyPressRequired, ClrScrRequired)
int KeyPressRequired, ClrScrRequired;
/*++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  List the contents of the normalisation array  */
/*++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	int i,j,k,fac,use;
	float ff,probval;
	char ch;
	float RowEffect[ResultHeight];

	if (ClrScrRequired == 1) ClearScreen();
	if (AnalysisLive == 0) 
          printf("Normalisation Array:- (LOADING)\n\n");
	if (AnalysisLive == 1) printf("Normalisation Array:- (LIVE)\n\n");
	for (i=0; i<ResultWidth; i++) printf(" %7.4f",Prob_Limit_Array[i]);
	printf("\n");
/*  Calculate Row multipliers  */
	fac = 1;
	for (i=0; i<ResultHeight; i++) {
		j = ResultHeight - 1 - i;
		fac = fac * 2;
		use = fac / 2;
		RowEffect[j] = use;
	}
/*  Do sum to build normalisation array  */
	for (i=0; i<ResultHeight; i++) {
		for (j=0; j<ResultWidth; j++) {
			probval = Prob_Limit_Array[j];
			ff = probval * NoRepeats * RowEffect[i];
			Normal[i][j] = ff;
			printf(" %7.1f",ff);
		}
		printf("\n");
	}
	printf("\r\n\n");
	if (KeyPressRequired == 1) ch = getch();
}



void ListScale(KeyPressRequired, ClrScrRequired)
int KeyPressRequired, ClrScrRequired;
/*+++++++++++++++++++++++++++++++++++++++++*/
/*  List the contents of the scaled array  */
/*+++++++++++++++++++++++++++++++++++++++++*/
{
	int i,j,k;
	float ff;
	char ch;

	if (ClrScrRequired == 1) ClearScreen();
	if (AnalysisLive == 0) 
          printf("Scale Array:- (LOADING)\n\n");
	if (AnalysisLive == 1) printf("Scale Array:- (LIVE)\n\n");
	for (i=0; i<ResultWidth; i++) 
          printf(" %7.4f",Prob_Limit_Array[i]);
	printf("\n");
/*  Do sum to build normalisation array  */
	for (i=0; i<ResultHeight; i++) {
		for (j=0; j<ResultWidth; j++) {
			ff = Result[i][j] / Normal[i][j];
			if (ff > 999.0) ff = 999.0;
			Scale[i][j] = ff;
			printf(" %7.3f",ff);
		}
		printf("\n");
	}
	printf("\r\n\n");
	if (KeyPressRequired == 1) ch = getch();
}



void ListAll()
/*+++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  List the contents of the analysis summary array  */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	char ch;

	ListAnalysisSummary(0,1);
	ListAnalysis(0,0);
	ListResults(0,0);
	ListNormalisation(0,0);
	ListScale(1,0);
}



void ClearResults()
/*+++++++++++++++++++++++++++++++*/
/*  Clear the Result data array  */
/*+++++++++++++++++++++++++++++++*/
{
	int i,j;

	for (i=0; i<ResultHeight; i++) {
		for (j=0; j<ResultWidth; j++) Result[i][j] = 0;
	}
}



void ClearAnalysis()
/*+++++++++++++++++++++++++++++++++*/
/*  Clear the Analysis data array  */
/*+++++++++++++++++++++++++++++++++*/
{
	int i,j;

	for (i=0; i<ResultHeight; i++) {
		for (j=0; j<AnalysisWidth; j++) Analysis[i][j] = 0;
	}
	for (j=0; j<ResultHeight; j++) {
		AnalysisMark[j] = 0;	/*  Counter for No of entries added  */
		AnalysisSum[j] = 0;	/*  Sum of entries added             */
		AnalysisBin[j] = 0;	/*  Address for next summation       */
		AnalysisCen[j] = HalfWidth;	/*  Burst analysis position  */
		AnalysisTry[j] = 0;	/*  No of burst bin trys	     */
	}
}



void FileWriteResults(num)
int num;
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Write the contents of the results array to a file         */
/*  Running set of result array tables so file open & close   */
/*  are not in this routine                                   */
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	int i,j,k,num1;

	num1 = num + 1;
	if (AnalysisLive == 0)
	fprintf(zzcon,"No. of Identified Bursts:- (LOADING)  Block %6d\n",num1);
  	if (AnalysisLive == 1) 
	fprintf(zzcon,"No. of Identified Bursts:- (LIVE)  Block %6d\n",num1);
	for (i=0; i<ResultWidth; i++) 
          fprintf(zzcon," %7.4f",Prob_Limit_Array[i]);
	fprintf(zzcon,"\n");
	for (i=0; i<ResultHeight; i++) {
		for (j=0; j<ResultWidth; j++) {
			k = Result[i][j];
			fprintf(zzcon," %7d",k);
		}
		fprintf(zzcon,"\n");
	}
	fprintf(zzcon,"\n");
}



void SaveLightCurve(num)
int num;
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Save a row from the Analysis Array - a light curve        */
/*  of AnalysisWidth bins                                     */
/*  The row used (SaveRowNo) is selected from the menu        */
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	int i,k,row;
        FILE * fp;

	fp = fopen("Saved_LC.dat","w");
	row = SaveRowNo - 1;
	fprintf(fp,"Saved Light Curve:-  ");
	fprintf(fp,"width = %5d, row = %3d, basic binsize = %5d\n",
	  AnalysisWidth,SaveRowNo,binsize);
	for (i=0; i<AnalysisWidth; i++) {
		k = Analysis[row][i];
		fprintf(fp,"%7d\n",k);
	}
	fclose(fp);
}