/* -*- C++ -*- */
/*************************************************************************
 * Copyright(c) 1995~2005  Masaharu Goto ([email protected])
 *
 * For the licensing terms see the file COPYING
 *
 ************************************************************************/
/*************************************************************************
* eular.c
*
* array pre-compiled class library is needed to run this demo program.
*************************************************************************/
#include <array.h>
const complex j=complex(0,1);
const double PI=3.141592;

main() {
  array x=array(-2*PI , 2*PI , 100 ); // start,stop,npoint
  plot << "Eular's Law" << x << exp(x*j) << "exp(x*j)" << endl ;
}

//.........这里部分代码省略.........
			if(ice_making <=ice_making_probability)
			{			
				check_icemaking = true;

				icemaker_running = true;

			//	ice_making_no = gl_random_sampled(3,ice_making_time);

				ice_making_no = 1;

				posted_power += icemaking_power;
				return_time = 60;			
				ice_making_no--;	

				if(ice_making_no == 0)
				{
					check_icemaking = false;
				}
			}
			else
			{
				icemaker_running = false;
				return_time = 60;	
			}
		}	
		else
		{
			icemaker_running = false;

		}
	}	
	
	load.power.SetPowerFactor(posted_power/1000, load.power_factor);
	load.admittance = complex(0,0,J); 
	load.current = complex(0,0,J);		

	if(true==door_open && true==door_to_open)
	{
		door_time = door_return_time<door_next_open_time?door_return_time:door_next_open_time;
	}
	else if(true==door_to_open)
	{
		door_time = door_next_open_time;
	}
	else if(true==door_open){
		
		door_time = door_return_time;

	}
	else if(start_time>0)
	{
		door_time = (3600 - ((t1-start_time)%3600));
	}

	if(0.0==return_time)
	{
		dt1 = cycle_time>door_time?door_time:cycle_time;		
	}
	else
	{
		if(cycle_time>return_time)
			if(return_time>door_time)
				if(door_time>long_compressor_cycle_time)
					dt1 = long_compressor_cycle_time;
				else
					dt1 = door_time;					

TIMESTAMP house::sync_panel(TIMESTAMP t0, TIMESTAMP t1)
{
	TIMESTAMP sync_time = TS_NEVER;
	OBJECT *obj = OBJECTHDR(this);

	// clear accumulators for panel currents
	complex I[3]; I[X12] = I[X23] = I[X13] = complex(0,0);

	// clear heatgain accumulator
	double heatgain = 0;

	// gather load power and compute current for each circuit
	CIRCUIT *c;
	for (c=panel.circuits; c!=NULL; c=c->next)
	{
		// get circuit type
		int n = (int)c->type;
		if (n<0 || n>2)
			GL_THROW("%s:%d circuit %d has an invalid circuit type (%d)", obj->oclass->name, obj->id, c->id, (int)c->type);
		/*	TROUBLESHOOT
			Invalid circuit types are an internal error for the house panel.  Please report this error.  The likely causes
			include an object that is not a house is being processed by the house model, or the panel was not correctly
			initialized.
		*/

		// if breaker is open and reclose time has arrived
		if (c->status==BRK_OPEN && t1>=c->reclose)
		{
			c->status = BRK_CLOSED;
			c->reclose = TS_NEVER;
			sync_time = t1; // must immediately reevaluate devices affected
			gl_debug("house:%d panel breaker %d closed", obj->id, c->id);
		}

		// if breaker is closed
		if (c->status==BRK_CLOSED)
		{
			// compute circuit current
			if ((c->pV)->Mag() == 0)
			{
				gl_debug("house:%d circuit %d (enduse %s) voltage is zero", obj->id, c->id, c->pLoad->name);
				break;
			}
			
			complex current = ~(c->pLoad->total*1000 / *(c->pV)); 

			// check breaker
			if (c->max_amps>0 && current.Mag()>c->max_amps)
			{
				// probability of breaker failure increases over time
				if (c->tripsleft>0 && gl_random_bernoulli(1/(c->tripsleft--))==0)
				{
					// breaker opens
					c->status = BRK_OPEN;

					// average five minutes before reclosing, exponentially distributed
					c->reclose = t1 + (TIMESTAMP)(gl_random_exponential(1/300.0)*TS_SECOND); 
					gl_debug("house:%d circuit breaker %d tripped - enduse %s overload at %.0f A", obj->id, c->id,
						c->pLoad->name, current.Mag());
				}

				// breaker fails from too frequent operation
				else
				{
					c->status = BRK_FAULT;
					c->reclose = TS_NEVER;
					gl_debug("house:%d circuit breaker %d failed", obj->id, c->id);
				}

				// must immediately reevaluate everything
				sync_time = t1; 
			}

			// add to panel current
			else
			{
				tload.power += c->pLoad->power;	// reminder: |a| + |b| != |a+b|
				tload.current += c->pLoad->current;
				tload.admittance += c->pLoad->admittance; // should this be additive? I don't buy t.a = c->pL->a ... -MH
				tload.total += c->pLoad->total;
				tload.heatgain += c->pLoad->heatgain;
				tload.energy += c->pLoad->power * gl_tohours(t1-t0);
				I[n] += current;
				c->reclose = TS_NEVER;
			}
		}

		// sync time
		if (sync_time > c->reclose)
			sync_time = c->reclose;
	}

	// compute line currents and post to meter
	if (obj->parent != NULL)
		LOCK_OBJECT(obj->parent);

	pLine_I[0] = I[X13];
	pLine_I[1] = I[X23];
	pLine_I[2] = 0;
	*pLine12 = I[X12];
//.........这里部分代码省略.........

bool analyzeTrafficBurst(signalVector &rxBurst,
			 unsigned TSC,
			 float detectThreshold,
			 int samplesPerSymbol,
			 complex *amplitude,
			 float *TOA,
			 unsigned maxTOA,
                         bool requestChannel,
                         signalVector **channelResponse,
			 float *channelResponseOffset) 
{

  assert(TSC<8);
  assert(amplitude);
  assert(TOA);
  assert(gMidambles[TSC]);

  if (maxTOA < 3*samplesPerSymbol) maxTOA = 3*samplesPerSymbol;
  unsigned spanTOA = maxTOA;
  if (spanTOA < 5*samplesPerSymbol) spanTOA = 5*samplesPerSymbol;

  unsigned startIx = (66-spanTOA)*samplesPerSymbol;
  unsigned endIx = (66+16+spanTOA)*samplesPerSymbol;
  unsigned windowLen = endIx - startIx;
  unsigned corrLen = 2*maxTOA+1;

  unsigned expectedTOAPeak = (unsigned) round(gMidambles[TSC]->TOA + (gMidambles[TSC]->sequenceReversedConjugated->size()-1)/2);

  signalVector burstSegment(rxBurst.begin(),startIx,windowLen);

  static complex staticData[200];
  signalVector correlatedBurst(staticData,0,corrLen);
  correlate(&burstSegment, gMidambles[TSC]->sequenceReversedConjugated,
					    &correlatedBurst, CUSTOM,true,
					    expectedTOAPeak-maxTOA,corrLen);

  float meanPower;
  *amplitude = peakDetect(correlatedBurst,TOA,&meanPower);
  float valleyPower = 0.0; //amplitude->norm2();
  complex *peakPtr = correlatedBurst.begin() + (int) rint(*TOA);

  // check for bogus results
  if ((*TOA < 0.0) || (*TOA > correlatedBurst.size())) {
        *amplitude = 0.0;
        return false;
  }

  int numRms = 0;
  for (int i = 2*samplesPerSymbol; i <= 5*samplesPerSymbol;i++) {
    if (peakPtr - i >= correlatedBurst.begin()) { 
      valleyPower += (peakPtr-i)->norm2();
      numRms++;
    }
    if (peakPtr + i < correlatedBurst.end()) {
      valleyPower += (peakPtr+i)->norm2();
      numRms++;
    }
  }

  if (numRms < 2) {
        // check for bogus results
        *amplitude = 0.0;
        return false;
  }

  float RMS = sqrtf(valleyPower/(float)numRms)+0.00001;
  float peakToMean = (amplitude->abs())/RMS;

  // NOTE: Because ideal TSC is 66 symbols into burst,
  //       the ideal TSC has an +/- 180 degree phase shift,
  //       due to the pi/4 frequency shift, that 
  //       needs to be accounted for.
  
  *amplitude = (*amplitude)/gMidambles[TSC]->gain;
  *TOA = (*TOA) - (maxTOA); 

  LOG(DEBUG) << "TCH peakAmpl=" << amplitude->abs() << " RMS=" << RMS << " peakToMean=" << peakToMean << " TOA=" << *TOA;

  LOG(DEBUG) << "autocorr: " << correlatedBurst;
  
  if (requestChannel && (peakToMean > detectThreshold)) {
    float TOAoffset = maxTOA; //gMidambles[TSC]->TOA+(66*samplesPerSymbol-startIx);
    delayVector(correlatedBurst,-(*TOA));
    // midamble only allows estimation of a 6-tap channel
    signalVector channelVector(6*samplesPerSymbol);
    float maxEnergy = -1.0;
    int maxI = -1;
    for (int i = 0; i < 7; i++) {
      if (TOAoffset+(i-5)*samplesPerSymbol + channelVector.size() > correlatedBurst.size()) continue;
      if (TOAoffset+(i-5)*samplesPerSymbol < 0) continue;
      correlatedBurst.segmentCopyTo(channelVector,(int) floor(TOAoffset+(i-5)*samplesPerSymbol),channelVector.size());
      float energy = vectorNorm2(channelVector);
      if (energy > 0.95*maxEnergy) {
	maxI = i;
	maxEnergy = energy;
      }
    }
	
    *channelResponse = new signalVector(channelVector.size());
    correlatedBurst.segmentCopyTo(**channelResponse,(int) floor(TOAoffset+(maxI-5)*samplesPerSymbol),(*channelResponse)->size());
//.........这里部分代码省略.........

      /*!
      \deprecated use standard contructors for typecasting

      \sa cxsc::complex::complex(const cdotprecision &)
      */
      friend inline complex _complex(const cdotprecision &a) throw() { return complex(a); }

scalar fn_init(double x, double y, scalar& dx, scalar& dy) {
  return complex(exp(-10*(x*x + y*y)), 0);
}

 /** exponentioal of a complex number
 \ingroup complex
 \param[in] z Complex number
 \return \f$ e^z \f$*/
 complex exp(const complex& z)
 {
   return complex(gsl_complex_exp(z.as_gsl_type()));
 }

 /** Logarithm of a complex number (base 10)
 \ingroup complex
 \param[in] z Complex number
 \return \f$ \log_{10} z \f$*/
 complex log10(const complex& z)
 {
   return complex(gsl_complex_log10(z.as_gsl_type()));
 }

 complex complex::conjugate() const
 {
     gsl_complex t = gsl_complex_conjugate(_complex);
     return complex(t);
 }

 complex complex::inverse() const
 {
     gsl_complex t = gsl_complex_inverse(_complex);
     return complex(t);
 }

 complex complex::operator/(const double& a) const
 {
   gsl_complex rl = gsl_complex_div_real(_complex,a);
   return complex(rl);
 }

 complex complex::operator/(const complex& z1) const
 {
   gsl_complex rl = gsl_complex_div(_complex, z1._complex);
   return complex(rl);
 }

 complex complex::operator-() const
   {
     gsl_complex t = gsl_complex_negative(_complex);
     return complex(&t);
   }

 /** Hyperbolic cosecant
 \ingroup complex
 \param[in] z Complex number
 \return \f$ \mathrm{csch} z \f$*/
 complex csch(const complex& z)
 {
   return complex(gsl_complex_csch(z.as_gsl_type()));
 }

 /** Logarithm of a complex number (base b)
 \ingroup complex
 \param[in] z Complex number
 \param[in] b Complex number
 \return \f$ \log_b z \f$*/
 complex log(const complex& z,
             const complex& b)
 {
   return complex(gsl_complex_log_b(z.as_gsl_type(),b.as_gsl_type()));
 }

 /** Inverse hyperbolic cotangent
 \ingroup complex
 \param[in] z Complex number
 \return \f$ \mathrm{acoth}(z) \f$*/
 complex arccoth(const complex& z)
 {
   return complex(gsl_complex_arccoth(z.as_gsl_type()));
 }

 /** DiLogarithm of a complex number
 \ingroup complex
 \param[in] z Complex number
 \return \f$ Li_2(z) \f$*/
 complex dilog(const complex& z)
 {
   gsl_sf_result re, im;
   gsl_sf_complex_dilog_xy_e(z.real(), z.imag(), &re, &im);
   return complex(re.val, im.val, false);
 }

 /** Square root of a complex number
 \ingroup complex
 \param[in] z Complex number
 \return \f$ \sqrt z \f$*/
 complex sqrt(const complex& z)
 {
   return complex(gsl_complex_sqrt(z.as_gsl_type()));
 }

complex bc_values(int marker, double x, double y)
{
  return complex(0.0, 0.0);
}

 /** Complex number to the z2 complex order
 \ingroup complex
 \param[in] z1 Complex number
 \param[in] z2 Complex number
 \return \f$ z_1^{z_2} \f$*/
 complex pow(const complex& z1,
             const complex& z2)
 {
   return complex(gsl_complex_pow(z1.as_gsl_type(),
                                  z2.as_gsl_type()));
 }