// quadratic interpolation: returns x of max, and the max itself in *r
double quint(double *r, double v1, double v2, double v3) ; 

void analysis(
			  LPFOURIER fou,	  // fourier struct
			        int N,		  // # of items in series
					int idx,	  // index of result item
			     double order,	  // order of required component
				    int Raw		  // which data series to analyse
			 )
{
      int i ;
   double * c, Ph, mPh, Tg ;
   double co, A, B, sum ;
   
     c = (double *)GlobalAlloc(GPTR, N * sizeof(double)) ;
   	Tg = align->period ;   // grand period

 // find the mean value of the data
if (Raw)
{	 
   for ( i = 0, sum = 0.0 ; i < N ; sum += fou->data[i++]) ;
        fou->DataMean = sum / (double)N ;
}
else
{
   for ( i = 0, sum = 0.0 ; i < N ; sum += fou->mdata[i++]) ;
        fou->mDataMean = sum / (double)N ;
}

 // establish the co-efficient for the
 // order of the required component

       co = twopi * order / (N-1) ;

 // The component A has a frequency (order) times higher
 // than that of the fundamental.  The period of the
 // fundamental is the grand period of the data series,
 // so the period of a component is (order) times smaller
 // than the grand period.
		
 // compute the cosine component amplitude
 // for the zero-mean data series at this order.

   for (i = 0 ; i < N ; i++)
   {
	   if (Raw)
     c[i] = (fou->data[i]  - fou->DataMean)  * cos(co * (double)i) ;
       else
     c[i] = (fou->mdata[i] - fou->mDataMean) * cos(co * (double)i) ;
   }
        A = 2.0 * integral(c, N) ;

 // compute the sine component amplitude
 // for the mean-extracted data series at this order
		
   for (i = 0 ; i < N ; i++)
   {
	   if (Raw)
     c[i] = (fou->data[i]  - fou->DataMean)  * sin(co * (double)i) ;
	   else
     c[i] = (fou->mdata[i] - fou->mDataMean) * sin(co * (double)i) ;
   }
        B = 2.0 * integral(c, N) ;

 // place results in the struct
   if (Raw)
   {
    fou->amp[idx]    = sqrt(A*A + B*B) ; // spectrum component amplitude
    fou->phase[idx]  = Ph = itan(B, A) ; // spectrum component phase
	fou->tShift[idx] = Tg * Ph /(twopi * order) ; // time shift
    fou->Ac[idx]     = A ;				//   Cosine component amplitude
    fou->Bs[idx]     = B ;				//     Sine component amplitude
   }
   else
   {
    fou->mamp[idx]    = sqrt(A*A + B*B) ;  // spectrum component amplitude
    fou->mphase[idx]  = mPh = itan(B, A) ; // spectrum component phase
	fou->mtShift[idx] = Tg * mPh /(twopi * order) ; // time shift
    fou->mAc[idx]     = A ;				   //   Cosine component amplitude
    fou->mBs[idx]     = B ;				   //     Sine component amplitude
   }

   GlobalFree(c) ;
}


void Spectrum( 
			  LPFOURIER fou,
			        int N,
			     double * pOrderStep,
			     double Hi_Order,
				 double Lo_Order,
				    int RawMean	  // 1 => Raw data, 0 => Mean data
			)
// N is the number of data points
// Lo_Order is the lowest order wanted.
// Hi_Order is the highest order wanted.
// OrderStep is the interval between successive orders wanted

// 4th June 2001 : OrderStep made a pointer so that any 
// adjustments made here propagate back to source. 
{
     int i ;
	 int mxori ;  // for test of MAXORDINDEX limit
  double order, ostep, temp, od ;
  static double f1 = 0.0, f2 = 0.0, f3 = 0.0 ;
  static double o1 = 0.0, o2 = 0.0, o3 = 0.0 ;
  static double p1 = 0.0, p2 = 0.0, p3 = 0.0 ;

         SendMessage(hWndPB1, PBM_SETPOS, (WPARAM)0, 0) ;

  // Adjust OrderStep to match
  // an integral number of
  // order-step indices.

  	temp = floor((od = Hi_Order - Lo_Order) /  * pOrderStep) ; 
	* pOrderStep = od / temp ;

	fou->OrderStep = ostep = * pOrderStep ;

  // OrderStep and index should now correspond exactly

    mxori = (int)ceil(od / ostep) ;

	if ( (mxori + 25) > MAXORDINDEX )
		MessageBox(NULL,"Number of Orders exceeds Limit",
		                "Too Many Fourier Orders",
						 MB_ICONEXCLAMATION|MB_OK
				  ) ;

	fou->max_order_index = mxori ; 

 // On this definition, indices are from zero to max_order_index inclusive.
 // So (max_order_index + 1) orders need to be computed.   

     for (
		             i = 0,
		   fou->Maxamp = fou->mMaxamp
		               = fou->bigPeak
					   = fou->mbigPeak
					   = -1.0e5 ;
	                i <= fou->max_order_index  ;
		            i++
		 )
      {

		order = (double)i * ostep + Lo_Order ;

        analysis( fou, N, i, order, RawMean) ;

		f1 = f2 ; o1 = o2 ; p1 = p2 ;
		f2 = f3 ; o2 = o3 ;	p2 = p3 ;
		f3 = (RawMean) ?   fou->amp[i] :   fou->mamp[i] ;
		p3 = (RawMean) ? fou->phase[i] : fou->mphase[i] ;
		o3 = order ;
		if (i > 2 && f1 < f2 && f2 > f3)      // we have a spectral peak
		{
			static double a, b, s[4] ;  // s[] used by quint()!!!

		  // the order of the peak	
            a = o1 + ostep * quint(s, f1, f2, f3) ;
			if (RawMean)
			    fou->ordPeak[fou->Peaks]   = a ;
			else
			    fou->mordPeak[fou->mPeaks] = a ;

		  // the value of the peak
			if (RawMean)
			    fou->spcPeak[fou->Peaks]   = s[0] ;
			else
			    fou->mspcPeak[fou->mPeaks] = s[0] ;

		  // the phase of the peak
			analysis( fou, N, MAXORDINDEX, a, RawMean ) ;
			b = (RawMean) ? fou->phase[MAXORDINDEX] : fou->mphase[MAXORDINDEX] ;
			if (RawMean)
				fou->phaPeak[fou->Peaks]   = b ;
			else
				fou->mphaPeak[fou->mPeaks] = b ;

		  // note properties of highest peak so far

			if (RawMean)
			{
			if (s[0] > fou->bigPeak) 
			    {
				 fou->bigPeak = s[0] ;
				 fou->bigOrd  = a ;
				 fou->bigPha  = b ;
			    }
			}
			else
			{
			if (s[0] > fou->mbigPeak)
				{
				 fou->mbigPeak = s[0] ;
				 fou->mbigOrd  = a ;
				 fou->mbigPha  = b ;
				}
			}

			if (RawMean)
			fou->Peaks++ ;
			else
			fou->mPeaks++ ;
		}

		if (RawMean)
		{
			if (fou->amp[i] >= fou->Maxamp)
			{
				fou->Maxamp = fou->amp[i] ;
				fou->Maxampidx = i ;
				fou->Maxampord = order ;
			}

		}
		else
		{
			if (fou->mamp[i] >= fou->mMaxamp)
			{
				fou->mMaxamp = fou->mamp[i] ;
				fou->mMaxampidx = i ;
				fou->mMaxampord = order ;
			}
		}

        SendMessage(hWndPB1,
		            PBM_SETPOS,
					(WPARAM)((UINT)(
					10000.0 * i / fou->max_order_index
					)), 0) ;

      }
			   if (RawMean) // is raw
			   {
				if (_is.mLamSpExists) // if we have the moving mean data as well, select greatest chart span
				fou->ComMaxamp = (fou->Maxamp > fou->mMaxamp) ? fou->Maxamp : fou->mMaxamp ;
				else				 // otherwise, use the current raw span
				fou->ComMaxamp = fou->Maxamp ;
			   }
			   else // is moving mean
			   {
				if (_is.LamSpExists) // if we have the raw data as well,  select greatest chart span
				fou->ComMaxamp = (fou->Maxamp > fou->mMaxamp) ? fou->Maxamp : fou->mMaxamp ;
				else				// otherwise, use the current moving mean span
				fou->ComMaxamp = fou->mMaxamp ;
			   }

        SendMessage(hWndPB1, PBM_SETPOS, (WPARAM)0, 0) ;
}

// Quadratic interpolation
double quint(double *r, double v1, double v2, double v3)
// Points (0,v1),(1,v2),(2,v3) assumed.
//
// Co-efficients a, b and c found for
// equation of parabola y = a*x*x + b*x + c.
//
// The y value at turning point at x,y returned in r[0]
// a, b and c are returned in r[1], r[2] and r[3] resp.
// The function returns the x value directly.
{
  double x, a, b, c ;
	r[1] = a = ( v3 + v1 )/2.0 - v2 ;
	r[2] = b = v2 - v1 - a ;
	   x =	-b / (2.0 * a) ;
	r[3] = c = v1 ;
    r[0] = a * x*x + b * x + c ;
	return x ;
}

void Synthesis(
			   LPFOURIER  f,		 // Data struct
			         int  N,		 // Number of values to synthesize
			      double  Hi_Order,	 // Highest Order of range
				  double  Lo_Order,	 // Lowest  Order of range
				  double   R_Order,	 // Resonant order
				  double       Sel,  // Selectivity factor
				  double      Gain,	 // Transfer Gain
					 int     Tuned,  // Whether tuned
					 int    Onbuds,	 // Whether on bud series or general dates
					 int       Raw,	 // Whether raw or mean data
					 int    Forced,  // Use forced co-efficients
					 int Alignment,	 // Synthesise Alignment data
					 int  Integral,  // Synthesize the integral
					 int Derivative  // Synthesize the derivative
			  )
{ 
      int i, j, k, * diptyp, * dp, * peaktyp, * pp ;
   double order, co, ostep, ospan ;
   double * c, * dipval, * peakval ;
   double h ;   
   double syn ;
   double Phase ;
   double ep1, period ;
   static double u, z, temp ;

   static double sr1, sr2, sr3 ;
   static double sm1, sm2, sm3 ;

   static double dr1, dr2, dr3 ;
   static double dm1, dm2, dm3 ;

   static dipcount, peakcount ;

	       k = f->max_order_index ;
           c = (double *)GlobalAlloc(GPTR, (k + 10) * sizeof(double)) ;

      dipval = (double *)GlobalAlloc(GPTR, (k + 10) * sizeof(double)) ;
          dp = (   int *)GlobalAlloc(GPTR, (k + 10) * sizeof(   int)) ;
	  diptyp = (   int *)GlobalAlloc(GPTR, (k + 10) * sizeof(   int)) ;

     peakval = (double *)GlobalAlloc(GPTR, (k + 10) * sizeof(double)) ;
          pp = (   int *)GlobalAlloc(GPTR, (k + 10) * sizeof(   int)) ;
	 peaktyp = (   int *)GlobalAlloc(GPTR, (k + 10) * sizeof(   int)) ;

	  ospan = Hi_Order - Lo_Order ;
	  ostep = ospan / (double)k ;
  diptyp[0] = 2 ;
	  dp[0] = 0 ;

  sr1 = sr2 = sr3 = sm1 = sm2 = sm3 = 0.0 ;
  dr1 = dr2 = dr3 = dm1 = dm2 = dm3 = 0.0 ;

  if (Tuned)
  {
	  f->maxtau = -10000.0 ;
	  f->mintau =  10000.0 ;

		  // this is for the selectivity graph
			     z = sqrt(3.0) / (2.0 * Sel) ;
				 u = sqrt(z*z + 1.0) ;
		 hi_dB_ord = fabs(R_Order / (z - u)) ;
		 lo_dB_ord = fabs(R_Order / (z + u)) ;
		 bandwidth = fabs(hi_dB_ord - lo_dB_ord) ;
		 if (hi_dB_ord < lo_dB_ord)
		 {
			  temp = hi_dB_ord ;
		 lo_dB_ord = hi_dB_ord ;
		 hi_dB_ord = temp ;
		 }
  }
  
  for(
	  i = 0,
	  dipcount = peakcount = 0 ; // set the counts to nil to stop useless summings
      i <= N ;
	  i++
	 )
   {
	double W, w, x;
		  
    for( j = 0 ; j <= k ; j++ )
     {
	  order = (double)j * ostep + Lo_Order  ;
         co = twopi * order/(double)(N-1) ;

		if (Tuned)
		{
				 if (order == 0.0)
				 w = 1e90 ;          // just make it huge
				 else
				 w = R_Order / order ;
			     W = 1.0 / w - w  ;
				 x = 1.0 / sqrt(1.0 + pow((Sel * W), 2.0)) ; // basic transfer function
			if (wantNotch)
			f->tau[j] = Gain * (1.0 - x) ;	// notch
			else
			f->tau[j] = Gain * x ;			// pass
		    f->eta[j] = atan(Sel * W) ;
			}// if (Tuned)

		if (Forced) // For now, just make all the amplitudes the same. 
		{
			  h = 0.025 ;
		}
		else
		{
		if (Raw)
		 h = sqrt(pow(f->Ac[j],2.0)  + pow(f->Bs[j],2.0)) ;
		else // running mean
		 h = sqrt(pow(f->mAc[j],2.0) + pow(f->mBs[j],2.0)) ;

		if (Tuned)
		  {
			h *= f->tau[j] ; // modify amplitudes by selectivity curve

			if (f->tau[j] > f->maxtau)
			{
				f->maxtau    = f->tau[j] ;
				f->taumaxord = order ;
				f->taumaxidx = j ;
				f->maxeta    = f->eta[j] ;
			}

			if (f->tau[j] < f->mintau)
			{
				f->mintau    = f->tau[j] ;
				f->tauminord = order ;
				f->tauminidx = j ;
				f->mineta    = f->eta[j] ;
			}
		  }// if (Tuned)
		}

		Phase = Raw ? f->phase[j] : f->mphase[j] ;

	    if (Tuned) Phase += f->eta[j] ; // modify phases by selectivity curve

		 c[j] = h * cos(co * (double)i - Phase) ;
     }

     if (Raw) // using raw lambdas
	 {
		 syn = ospan * integral(c, k) + f->DataMean ;  // synthesised value

	  if (wantBumps)  // looking for peaks or dips in the synthesis
	  {
	  sr1 = sr2 ; sr2 = sr3 ; sr3 = syn ;  // left shift old and assign new  
	  dr1 = dr2 ; dr2 = dr3 ; dr3 = syn ;

// DIPS

	  if (i && (!_is.RawDipsExist || (Tuned && !_is.RawTunedDipsExist))) // either tuned or untuned
	  {
	  if (sr1 > sr2 && sr2 == sr3) // falling to plateau
	  {
		  diptyp[dipcount  ] = 1 ;
		  dipval[dipcount  ] = sr2 ;
		      dp[dipcount++] = i - 1 ;
	  }
	  if (sr1 == sr2 && sr2 < sr3) // rising from plateau
	  {
		  diptyp[dipcount  ] = 2 ;
		  dipval[dipcount  ] = sr2 ;
		      dp[dipcount++] = i - 1 ;
	  }
	  if (sr1 > sr2 && sr2 < sr3) // single dip
	  {
		  diptyp[dipcount  ] = 3 ;
		  dipval[dipcount  ] = sr2 ;
		      dp[dipcount++] = i - 1 ;
	  }
	  } // end if (some raw dips don't exist)

// PEAKS

	  if (i && (!_is.RawPeaksExist || (Tuned && !_is.RawTunedPeaksExist))) // either tuned or untuned
	  {
	  if (dr1 < dr2 && dr2 == dr3) // rising to plateau
	  {
		  diptyp[peakcount  ] = 1 ;
		  dipval[peakcount  ] = dr2 ;
		      pp[peakcount++] = i - 1 ;
	  }
	  if (dr1 == dr2 && dr2 > dr3) // falling from plateau
	  {
		  peaktyp[peakcount  ] = 2 ;
		  peakval[peakcount  ] = dr2 ;
		       pp[peakcount++] = i - 1 ;
	  }
	  if (dr1 < dr2 && dr2 > dr3) // single peak
	  {
		  peaktyp[peakcount  ] = 3 ;
		  peakval[peakcount  ] = dr2 ;
		       pp[peakcount++] = i - 1 ;
	  }
	  } // if (some raw peaks don't exist)
	  } // if (wantBumps)

 	  if (Forced)
	  {
		  if (Alignment)
		  f->fsaldata[i] = syn ; // forced alignment
	      else
		    f->fsdata[i] = syn ; // forced raw lambda
	  }
	  else // standard, unforced synthesis
	  {
		  if (Tuned)
		  f->tsdata[i] = syn ;  // tuned raw lambda
		  else // untuned
		  {
			  if (wantOneDay)
				f->sOneDay[i] = syn ; // untuned daily raw lambda
			    else
				f->sdata[i] =  syn ;  // untuned raw lambda
		  }
	  }
	 }
	 else	// using running mean lambdas
	 {
		 syn = ospan * integral(c, k) + f->mDataMean ; // synthesised value
	 
	  if (wantBumps)
	  {
	  sm1 = sm2 ; sm2 = sm3 ; sm3 = syn ;  // left shift old and assign new
	  dm1 = dm2 ; dm2 = dm3 ; dm3 = syn ;

// DIPS

	  if (i &&(!_is.MeanDipsExist || (Tuned && !_is.MeanTunedDipsExist)))
	  {
	  if (sm1 > sm2 && sm2 == sm3) // type 1, falling to plateau
	  {
		  diptyp[dipcount  ] = 1 ;
		  dipval[dipcount  ] = sm2 ;
		      dp[dipcount++] = i - 1 ;
	  }
	  if (sm1 == sm2 && sm2 < sm3) // type 2, rising from plateau
	  {
		  diptyp[dipcount  ] = 2 ;
		  dipval[dipcount  ] = sm2 ;
		      dp[dipcount++] = i - 1 ;
	  }
	  if (sm1 > sm2 && sm2 < sm3)  // type 3, single dip
	  {
		  diptyp[dipcount  ] = 3 ;
		  dipval[dipcount  ] = sm2 ;
		      dp[dipcount++] = i - 1 ;
	  }
	  } // end if (some mean dips don't exist)

// PEAKS

	  if (i && (!_is.MeanPeaksExist || (Tuned && !_is.MeanTunedPeaksExist))) // either tuned or untuned
	  {
	  if (dm1 < dm2 && dm2 == dm3) // rising to plateau
	  {
		  diptyp[peakcount  ] = 1 ;
		  dipval[peakcount  ] = dm2 ;
		      pp[peakcount++] = i - 1 ;
	  }
	  if (dm1 == dm2 && dm2 > dm3) // falling from plateau
	  {
		  peaktyp[peakcount  ] = 2 ;
		  peakval[peakcount  ] = dm2 ;
		       pp[peakcount++] = i - 1 ;
	  }
	  if (dm1 < dm2 && dm2 > dm3) // single peak
	  {
		  peaktyp[peakcount  ] = 3 ;
		  peakval[peakcount  ] = dm2 ;
		       pp[peakcount++] = i - 1 ;
	  }
	  } // end if (some mean peaks don't exist)
	  } // end if (wantBumps)

	  if (Forced)
	  {
		  if (Alignment)
		  f->fsaldata[i] = syn ;	// same as under Raw, above: forced alignment
	      else
		   f->fmsdata[i] = syn ;	// forced mean lambda
	  }
	  else // standard, unforced synthesis
	  {
		  if (Tuned)
		  f->tmsdata[i] = syn ;	// tuned mean lambda
		  else // untuned
		  {
			  if (wantOneDay)
				  f->msOneDay[i] = syn ;	// untuned daily mean lambda
				  else
				  f->msdata[i]   = syn ;	// untuned mean lambda
		  }
	  }
	 } // using running mean lambda
   }// end of major loop

   if (wantBumps)
   {
   // sort out and record dips and peaks

	if (Onbuds)
	{
	   ep1 = BudEp1 ; // on bud series dates
	period = BudGP ;
	}
	else
	{
	   ep1 = AstEp1 ; // on general dates
	period = AstGP ;
	}

  // first clear out any old dip or peak data from the structs.
  // the negative flags signify that we are about to (re)create the data
  // and want to be sure of a clear field
   if (Raw) 
   {
	   if (!_is.RawDipsExist) // clear old raw dip data
	   {
	   for (i = 0 ; i < MAXPEAKS ; i++)
	   f->DipEps[i] = f->DipVal[i] = 0.0 ;
	   f->nDips  = 0 ;
	   }

	   if (!_is.RawPeaksExist) // clear old raw peak data
	   {
	   for (i = 0 ; i < MAXPEAKS ; i++)
	   f->PeakEps[i] = f->PeakVal[i] = 0.0 ;
	   f->nPeaks  = 0 ;
	   }

	if (Tuned && !_is.RawTunedDipsExist) // clear raw tuned dip data
	{
	   for (i = 0 ; i < MAXPEAKS ; i++)
	   f->tDipEps[i] = f->tDipVal[i] = 0.0 ;
	   f->ntDips = 0 ;
	}

	if (Tuned && !_is.RawTunedPeaksExist) // clear raw tuned peak data
	{
	   for (i = 0 ; i < MAXPEAKS ; i++)
	   f->tPeakEps[i] = f->tPeakVal[i] = 0.0 ;
	   f->ntPeaks = 0 ;
	}
   }
   else // is on mean data
   {
	   if (!_is.MeanDipsExist) // clear out old mean dip data
	   {
	   for (i = 0 ; i < MAXPEAKS ; i++)
	   f->mDipEps[i] = f->mDipVal[i] = 0.0 ;
	   f->nmDips  = 0 ;
	   }

	   if (!_is.MeanPeaksExist) // clear out old mean peak data
	   {
	   for (i = 0 ; i < MAXPEAKS ; i++)
	   f->mPeakEps[i] = f->mPeakVal[i] = 0.0 ;
	   f->nmPeaks  = 0 ;
	   }

    if (Tuned && !_is.MeanTunedDipsExist) // initialise or clear out old mean tuned dip data
	{
	   for (i = 0 ; i < MAXPEAKS ; i++)
	   f->tmDipEps[i] = f->tmDipVal[i] = 0.0 ;
	   f->nmtDips = 0 ;
	}

    if (Tuned && !_is.MeanTunedPeaksExist) // initialise or clear out old mean tuned peak data
	{
	   for (i = 0 ; i < MAXPEAKS ; i++)
	   f->tmPeakEps[i] = f->tmPeakVal[i] = 0.0 ;
	   f->nmtPeaks = 0 ;
	}
   }

// DIPS

   if (dipcount) // if we have any dips at all
   {
   for ( i = 0 ; i <= (dipcount + 1) ; i++ )
   {
    if (diptyp[i] == 2) // if current is rising from plateau...
    {
     if (diptyp[i - 1] == 1) // then if previous was falling to plateau....
	 {
		 int d ;
		// and there is one interval, or two intervals, between them, ....
		 if (
			 ((d = dp[i] - dp[i - 1]) == 1) || (d == 2)
			)
		{

		// ...then we have either a pair of adjacent "low corners", or
		// a pair of near-adjacent "low corners", so split the difference.

		if (Raw)
		{
			if (Tuned)
			{
			f->tDipVal[f->ntDips] = dipval[i-1] ;
			f->tDipEps[f->ntDips++] =
				ep1 + period * ((double)(dp[i-1] + dp[i]) / 2.0) / (double)N ;
			}
			else
			{
			f->DipVal[f->nDips] = dipval[i-1] ;
			f->DipEps[f->nDips++] =
				ep1 + period * ((double)(dp[i-1] + dp[i]) / 2.0) / (double)N ;
			}
		}
		else // is running mean
		{
			if (Tuned)
			{
			f->tmDipVal[f->nmtDips] = dipval[i-1] ;
			f->tmDipEps[f->nmtDips++] =
				ep1 + period * ((double)(dp[i-1] + dp[i]) / 2.0) / (double)N ;
			}
			else
			{
			f->mDipVal[f->nmDips] = dipval[i-1] ;
			f->mDipEps[f->nmDips++] =
				ep1 + period * ((double)(dp[i-1] + dp[i]) / 2.0) / (double)N ;
			}
		}
	    }
	 }
	}

// However,
   if (diptyp[i] == 3)
	{
   // then we have a "simple dip". Use it directly.
	   if (Raw)
	   {
		   if(Tuned)
		   {
			f->tDipVal[f->ntDips] = dipval[i] ;
			f->tDipEps[f->ntDips++] =
				ep1 + period * (double)dp[i] / (double)N ;
		   }
		   else
		   {
			f->DipVal[f->nDips] = dipval[i] ;
			f->DipEps[f->nDips++] =
				ep1 + period * (double)dp[i] / (double)N ;
		   }
	   }
	   else
	   {
		   if(Tuned)
		   {
			f->tmDipVal[f->nmtDips] = dipval[i] ;
 			f->tmDipEps[f->nmtDips++] =
				ep1 + period * (double)dp[i] / (double)N ;
		   }
		   else
		   {
			f->mDipVal[f->nmDips] = dipval[i] ;
 			f->mDipEps[f->nmDips++] =
				ep1 + period * (double)dp[i] / (double)N ;
		   }
	   }
   }
   }

   if ( i > dipcount + 1) // have we left the loop?
   {
   if (Raw && !Tuned)
   {
	   static n, count[20], idx, topidx = -1, botidx = 100 ;
	   double t, d, D ;
	   for (i = 0 ; i < 20 ; count[i++] = 0) ;

				for (i = 0 ; i <= f->nDips ; i++)
				{
					for (n = 0 ; n <= align->found ; n++)
					{
						if ((t = Date_Time[UsePlanet][n]->epal) >= (d = f->DipEps[i]))
						{
							f->rdpShift[i] = (d - t) ;	// get the dip shift from ordinary alignment
							break ;
						}
					}
				}

				// build simple Agraph, counting backwards from alignment epochs
				// to limit I.A.P. dependency
					for ( i = 0 ; i < f->nDips ; i++ )
					{
						idx = (int)fabs(f->rdpShift[i]) ;			// dip shift to nearest day as index
						count[idx]++ ;								// tot up dips falling in this day
						topidx = (idx >= topidx) ? idx : topidx ;	// earliest day in any IAP
						botidx = (idx <= botidx) ? idx : botidx ;	//   latest day in any IAP
					}

				// find and count temporal correlation clusters


				// find the temporal centroid of the series
					for ( i = 0, D = 0.0 ; i < f->nDips ; i++ ) D += f->DipEps[i] ;		// sum all the dip epochs
					f->rdpMeanEpoch = D / f->nDips ;		// this is the mean epoch

					_is.RawDipsExist = 1 ;
	   EnableWindow(GetDlgItem(hDlgDisplay,    IDC_C_RAWDIPS),	 TRUE) ;
	   EnableWindow(GetDlgItem(hRastDlg, IDC_C_SYNTH_RAW_DIPS),  TRUE) ;
   }

   if (Raw && Tuned)
   {
	   int n ;
	   double t, d ;
				for (i = 0 ; i <= f->ntDips ; i++)
				{
					for (n = 0 ; n <= align->found ; n++)
					{
						if ((t = Date_Time[UsePlanet][n]->epal) >= (d = f->tDipEps[i]))
						{
							f->trdpShift[i] = d - t ; break ;
						}
					}
				}
		_is.RawTunedDipsExist = 1 ;
		EnableWindow(GetDlgItem(hDlgDisplay, IDC_C_RAWTUNEDDIPS),    TRUE) ;
	    EnableWindow(GetDlgItem(hRastDlg,	 IDC_C_TUNED_RAW_DIPS),  TRUE) ;
   }

   if (!Raw && !Tuned)
   {
	   int n ;
	   double t, d ;
				for (i = 0 ; i <= f->nmDips ; i++)
				{
					for (n = 0 ; n <= align->found ; n++)
					{
						if ((t = Date_Time[UsePlanet][n]->epal) >= (d = f->mDipEps[i]))
						{
							f->mdpShift[i] = d - t ; break ;
						}
					}
				}
		_is.MeanDipsExist = 1 ;
		EnableWindow(GetDlgItem(hDlgDisplay,    IDC_C_MEANDIPS),         TRUE) ;
	    EnableWindow(GetDlgItem(hRastDlg,		IDC_C_SYNTH_MEAN_DIPS),  TRUE) ;
   }

   if (!Raw && Tuned)
   {
	   int n ;
	   double t, d ;
				for (i = 0 ; i <= f->nmtDips ; i++)
				{
					for (n = 0 ; n <= align->found ; n++)
					{
						if ((t = Date_Time[UsePlanet][n]->epal) >= (d = f->tmDipEps[i]))
						{
							f->tmdpShift[i] = d - t ; break ;
						}
					}
				}
		_is.MeanTunedDipsExist = 1 ;
		EnableWindow(GetDlgItem(hDlgDisplay,    IDC_C_MEANTUNEDDIPS),    TRUE) ;
	    EnableWindow(GetDlgItem(hRastDlg,		IDC_C_TUNED_MEAN_DIPS),  TRUE) ;
   }
   }

   } // if (dipcount)

   // PEAKS

   if (peakcount) // if we have any peaks at all
   {
   for ( i = 0 ; i <= (peakcount + 1) ; i++ )
   {
    if (peaktyp[i] == 2) // if current is falling from plateau
    {
     if (peaktyp[i - 1] == 1) // then if previous was falling to plateau
	 {
		 int d ;
		// and there is one interval, or two intervals, between them,
		 if (
			 ((d = pp[i] - pp[i - 1]) == 1)
			    ||
			 (d == 2)
			)
		{

		// then we have either a pair of adjacent "high corners", or
		// a pair of near-adjacent "high corners", so split the difference.

		if (Raw)
		{
			if (Tuned)
			{
			f->tPeakVal[f->ntPeaks] = peakval[i-1] ;
			f->tPeakEps[f->ntPeaks++] =
				ep1 + period * ((double)(pp[i-1] + pp[i]) / 2.0) / (double)N ;
			}
			else
			{
			f->PeakVal[f->nPeaks] = peakval[i-1] ;
			f->PeakEps[f->nPeaks++] =
				ep1 + period * ((double)(pp[i-1] + pp[i]) / 2.0) / (double)N ;
			}
		}
		else // is running mean
		{
			if (Tuned)
			{
			f->tmPeakVal[f->nmtPeaks] = peakval[i-1] ;
			f->tmPeakEps[f->nmtPeaks++] =
				ep1 + period * ((double)(pp[i-1] + pp[i]) / 2.0) / (double)N ;
			}
			else
			{
			f->mPeakVal[f->nmPeaks] = peakval[i-1] ;
			f->mPeakEps[f->nmPeaks++] =
				ep1 + period * ((double)(pp[i-1] + pp[i]) / 2.0) / (double)N ;
			}
		}
	    }
	 }
	}

// However,
   if (peaktyp[i] == 3)
	{
   // then we have a "simple peak". Use it directly.
	   if (Raw)
	   {
		   if(Tuned)
		   {
			f->tPeakVal[f->ntPeaks] = peakval[i] ;
			f->tPeakEps[f->ntPeaks++] =
				ep1 + period * (double)pp[i] / (double)N ;
		   }
		   else
		   {
			f->PeakVal[f->nPeaks] = peakval[i] ;
			f->PeakEps[f->nPeaks++] =
				ep1 + period * (double)pp[i] / (double)N ;
		   }
	   }
	   else
	   {
		   if(Tuned)
		   {
			f->tmPeakVal[f->nmtPeaks] = peakval[i] ;
 			f->tmPeakEps[f->nmtPeaks++] =
				ep1 + period * (double)pp[i] / (double)N ;
		   }
		   else
		   {
			f->mPeakVal[f->nmPeaks] = peakval[i] ;
 			f->mPeakEps[f->nmPeaks++] =
				ep1 + period * (double)pp[i] / (double)N ;
		   }
	   }
   }
   }

   if ( i > peakcount + 1) // have we left the loop?
   {
   if (Raw && !Tuned)
   {
	   int n ;
	   double t, d ;
				for (i = 0 ; i <= f->nPeaks ; i++)
				{
					for (n = 0 ; n <= align->found ; n++)
					{
						if ((t = Date_Time[UsePlanet][n]->epal) >= (d = f->PeakEps[i]))
						{
							f->rpkShift[i] = d - t ; break ;
						}
					}
				}
	   _is.RawPeaksExist = 1 ;
	   EnableWindow(GetDlgItem(hDlgDisplay, IDC_C_RAWPEAKS),        TRUE) ;
	   EnableWindow(GetDlgItem(hRastDlg,	IDC_C_SYNTH_RAW_PEAKS), TRUE) ;
   }

   if (Raw && Tuned)
   {
	   int n ;
	   double t, d ;
				for (i = 0 ; i <= f->ntPeaks ; i++)
				{
					for (n = 0 ; n <= align->found ; n++)
					{
						if ((t = Date_Time[UsePlanet][n]->epal) >= (d = f->tPeakEps[i]))
						{
							f->trpkShift[i] = d - t ; break ;
						}
					}
				}
		_is.RawTunedPeaksExist = 1 ;
		EnableWindow(GetDlgItem(hDlgDisplay,    IDC_C_RAWTUNEDPEAKS),   TRUE) ;
	    EnableWindow(GetDlgItem(hRastDlg,		IDC_C_TUNED_RAW_PEAKS), TRUE) ;
   }

   if (!Raw && !Tuned)
   {
	   int n ;
	   double t, d ;
				for (i = 0 ; i <= f->nmPeaks ; i++)
				{
					for (n = 0 ; n <= align->found ; n++)
					{
						if ((t = Date_Time[UsePlanet][n]->epal) >= (d = f->mPeakEps[i]))
						{
							f->mpkShift[i] = d - t ; break ;
						}
					}
				}
	   _is.MeanPeaksExist = 1 ;
	   EnableWindow(GetDlgItem(hDlgDisplay, IDC_C_MEANPEAKS),        TRUE) ;
	   EnableWindow(GetDlgItem(hRastDlg,	IDC_C_SYNTH_MEAN_PEAKS), TRUE) ;
   }

   if (!Raw && Tuned)
   {	   
	   int n ;
	   double t, d ;
				for (i = 0 ; i <= f->nmtPeaks ; i++)
				{
					for (n = 0 ; n <= align->found ; n++)
					{
						if ((t = Date_Time[UsePlanet][n]->epal) >= (d = f->tmPeakEps[i]))
						{
							f->tmpkShift[i] = d - t ; break ;
						}
					}
				}
	   _is.MeanTunedPeaksExist = 1 ;
	   EnableWindow(GetDlgItem(hDlgDisplay, IDC_C_MEANTUNEDPEAKS),   TRUE) ;
	   EnableWindow(GetDlgItem(hRastDlg,	IDC_C_TUNED_MEAN_PEAKS), TRUE) ;
   }
   }

   } // if (peakcount)
   } // if (wantBumps)

// finished with this lot,
// release the storage.
   GlobalFree(c) ;
   GlobalFree(dipval) ;
   GlobalFree(dp) ;
   GlobalFree(diptyp) ;
   GlobalFree(peakval) ;
   GlobalFree(pp) ;
   GlobalFree(peaktyp) ;
    
}