void BoatDetector::houghline(IplImage* edgeImage, IplImage* image, IplImage* lineImage) {

	//validation
	int points = 50; // points per row
	int rows = 3; // number of rows
	int ver_dist = 10; // vertical distance between points
	int hor_dist = image->width / points; // horizontal distance between points

	cvCopy(edgeImage, lineImage);

	CvSeq* hough_lines = 0;

	CvScalar line_color = cvScalar(120);

	hough_lines = cvHoughLines2( edgeImage, hough_storage, CV_HOUGH_STANDARD, 1, CV_PI/180, 100, 0, 0 );


	if(hough_lines->total == 0) {
		return;
	}

	bool find = false;

	CvPoint pt1, pt2;
	float* line;
	float theta;
	float rho;
	double a, b, x0, y0;
	for( int i = 0; i < min(hough_lines->total, 100); i++ )
	{
		line = (float*)cvGetSeqElem(hough_lines, i);
		theta = line[1];
		if(theta < 1.50 || theta > 1.60) {
			continue;
		}
		rho = line[0];

		a = cos(theta);
		b = sin(theta);
		x0 = a*rho;
		y0 = b*rho;
		pt1.x = cvRound(x0 + 1000*(-b));
		pt1.y = cvRound(y0 + 1000*(a));
		pt2.x = cvRound(x0 - 1000*(-b));
		pt2.y = cvRound(y0 - 1000*(a));		

		cvLine( lineImage, pt1, pt2, line_color, 2, CV_AA, 0 );
		find = true;
	}
	if(!find) {
		return;
	}

	bool run = true;
	int search_limit = lineImage->height - (ver_dist * rows);
	int line_step = lineImage->widthStep/sizeof(char);
	int img_step = image->widthStep/sizeof(char);
	int max_left, max_right;
	int tmp_limit;
	double count;
	while(run) {
		max_left = 0;
		max_right = 0;
		for(int i = ver_dist * rows; i < search_limit; i++) {
			if(((uchar)lineImage->imageData[i*line_step+3]) == line_color.val[0]) {
				if(i > max_left) {
					max_left = i;
				}
			}
			if(((uchar)lineImage->imageData[i*line_step + (lineImage->width-3)]) == line_color.val[0]) {
				if(i > max_right) {
					max_right = i;
				}
			}		
		}
		if(max_left == 0 || max_right == 0) {
			run = false;
			continue;
		}

		tmp_limit = (max_left + max_right) / 2;

		//limit validation
		count = 0;
		if(abs(max_left - max_right) < 10) {

			for(int i = tmp_limit - (ver_dist * rows), k = 0, t = rows*2; k < rows; i+=ver_dist, k++, t-=2) {
				for(int j = hor_dist; j < image->width; j+=hor_dist) {
					if(abs(image->imageData[i*img_step + j] -
						image->imageData[(i+t*ver_dist)*img_step + j] ) > 10 )
					{
						count++;
					}
				}
			}
		}
		if((count / (points * rows)) > 0.9 ) {

			sea_limit = tmp_limit;

//.........这里部分代码省略.........

//.........这里部分代码省略.........
                   exit (3);
                   }
                imag=(double) ntohs(data_short);
               
             }
             else
             {
                if (fread(&data_int,no_bytes,1,in_file) != 1)
                   {
                   fprintf (stderr, "Error in reading input data (int/float)\n");
                   exit (3);
                   }
                uval.ival = ntohl(data_int);
                if (bit[3] == 0)
                   real = (double) uval.ival;
                else
                   real = (double) uval.fval;
                if (fread(&data_int,no_bytes,1,in_file) != 1)
                   {
                   fprintf (stderr, "Error in reading input data (int/float)\n");
                   exit (3);
                   }
                uval.ival = ntohl(data_int);
                if (bit[3] == 0)
                   imag = (double) uval.ival;
                else
                   imag = (double) uval.fval;
             }             

             phi=atan2(imag,real);

             phi=phi+2.0* M_PI *delta_f*(j-lsfid)/sw;

             abs_data=sqrt(real*real+imag*imag);

             real=abs_data*cos(phi);  
             imag=abs_data*sin(phi);   

             if (no_bytes == 2)
                {
                data_short = (short)real;
                data_short = htons(data_short);
                if (fwrite(&data_short,no_bytes,1,out_file) != 1)
                   {
                   fprintf (stderr, "Error in writting output data (short)\n");
                   exit (3);
                   }     
 
                data_short = (short)imag;
                data_short = htons(data_short);
                if (fwrite(&data_short,no_bytes,1,out_file) != 1)
                   {
                   fprintf (stderr, "Error in writting output data (short)\n");
                   exit (3);
                   }
                }
              else
                {
                     
                if (bit[3] == 0)
                   uval.ival = (int)real;
                else
                   uval.fval = (float) real;
                data_int = htonl(uval.ival);
                if (fwrite(&data_int,no_bytes,1,out_file) != 1)
                   {
                   fprintf (stderr, "Error in writting output data (int/float)\n");
                   exit (3);
                   }     
                if (bit[3] == 0)
                   uval.ival = (int)imag;
                else
                   uval.fval = (float) imag;
                data_int = htonl(uval.ival);
 
                if (fwrite(&data_int,no_bytes,1,out_file) != 1)
                   {
                   fprintf (stderr, "Error in writting output data (int/float)\n");
                   exit (3);
                   }          
                }
            }
        }
   fprintf (stderr, "\nB0 CORRECTION WAS PERFORMED!\n");
   
/*** CLOSING FILES *************************************/                                 

      if (fclose(in_file) != 0) 
         {
         printf ("Error closing input file"); 
         exit (3);      
         }

      if (fclose(out_file) != 0) 
         {
         printf ("Error closing output file"); 
         exit (3);      
         }   
   exit(EXIT_SUCCESS);
}

void
set_M(QLA_ColorMatrix *m, int i)
{
#if 0
  static QLA_ColorMatrix t;
  for(int j=0; j<QLA_Nc; j++) {
    for(int k=0; k<QLA_Nc; k++) {
      QLA_c_eq_r_plus_ir(QLA_elem_M(*m,j,k),
			 (((j-k+QLA_Nc+1)*(j+k+1))%19)+cos(i),
			 (((j+4)*(k+1))%17)+sin(i));
      //QLA_real(QLA_elem_M(*m,j,k)) = 1;
      //QLA_imag(QLA_elem_M(*m,j,k)) = 0;
    }
  }
#endif
  for(int j=0; j<QLA_Nc; j++) {
    for(int k=0; k<QLA_Nc; k++) {
      QLA_c_eq_r(QLA_elem_M(*m,j,k), 0);
    }
  }
  QLA_Real step = 1e-5;
  if(Mtype&MtypeNZ) {
    for(int j=0; j<QLA_Nc; j++) {
      QLA_c_peq_r_plus_ir(QLA_elem_M(*m,j,j), step, -step);
    }
  }
  int ii=i;
  if((Mtype&MtypeNN)==0) ii>>=QLA_Nc;
  for(int j=0,k=1; ii; ii>>=1,j++) {
    if(j>=QLA_Nc) { j=0; k*=2; }
    if(ii&1) QLA_c_peq_r_plus_ir(QLA_elem_M(*m,j,j), k*step, -k*step);
  }
  ii = i;
  if((Mtype&MtypeNN)==0) {
    for(int j=0; j<QLA_Nc; j++) {
      if(ii&1) QLA_c_eqm_c(QLA_elem_M(*m,j,j), QLA_elem_M(*m,j,j));
      ii >>= 1;
    }
  }
  if(Mtype&MtypeH) { // make Hermitian
    QLA_ColorMatrix m2;
    QLA_M_eq_M(&m2, m);
    QLA_M_peq_Ma(&m2, m);
    QLA_M_eq_M(m, &m2);
  }
  if((Mtype&MtypeP)&&(Mtype&MtypeH)) { // make positive Hermitian
    QLA_ColorMatrix m2;
    QLA_M_eq_M_times_Ma(&m2, m, m);
    QLA_M_eq_M(m, &m2);
  }
  if(Mtype&MtypeA) { // make anti-Hermitian
    QLA_ColorMatrix m2;
    QLA_M_eq_M(&m2, m);
    QLA_M_meq_Ma(&m2, m);
    QLA_M_eq_M(m, &m2);
  }
  if((Mtype&MtypeT)&&(Mtype&MtypeA)) { // make traceless anti-Hermitian
    QLA_ColorMatrix m2;
    QLA_M_eq_antiherm_M(&m2, m);
    QLA_M_eq_M(m, &m2);
  }
  //QLA_Real n2;
  //QLA_r_eq_norm2_M(&n2, m);
  //printf("%i\t%g\n", i, n2);
}

static void HandleEntityAdjustment(void)
{
	char		*value;
	vec3_t		origin, newOrigin, angles;
	char		temp[MAX_QPATH];
	float		rotation;

	G_SpawnString("origin", NOVALUE, &value);
	if (Q_stricmp(value, NOVALUE) != 0)
	{
		if ( sscanf( value, "%f %f %f", &origin[0], &origin[1], &origin[2] ) != 3 ) {
			G_Printf( "HandleEntityAdjustment: failed sscanf on 'origin' (%s)\n", value );
			VectorClear( origin );
		}
	}
	else
	{
		origin[0] = origin[1] = origin[2] = 0.0;
	}

	rotation = DEG2RAD(level.mRotationAdjust);
	newOrigin[0] = origin[0]*cos(rotation) - origin[1]*sin(rotation);
	newOrigin[1] = origin[0]*sin(rotation) + origin[1]*cos(rotation);
	newOrigin[2] = origin[2];
	VectorAdd(newOrigin, level.mOriginAdjust, newOrigin);
	// damn VMs don't handle outputing a float that is compatible with sscanf in all cases
	Com_sprintf(temp, sizeof( temp ), "%0.0f %0.0f %0.0f", newOrigin[0], newOrigin[1], newOrigin[2]);
	AddSpawnField("origin", temp);

	G_SpawnString("angles", NOVALUE, &value);
	if (Q_stricmp(value, NOVALUE) != 0)
	{
		if ( sscanf( value, "%f %f %f", &angles[0], &angles[1], &angles[2] ) != 3 ) {
			G_Printf( "HandleEntityAdjustment: failed sscanf on 'angles' (%s)\n", value );
			VectorClear( angles );
		}

		angles[YAW] = fmod(angles[YAW] + level.mRotationAdjust, 360.0f);
		// damn VMs don't handle outputing a float that is compatible with sscanf in all cases
		Com_sprintf(temp, sizeof( temp ), "%0.0f %0.0f %0.0f", angles[0], angles[1], angles[2]);
		AddSpawnField("angles", temp);
	}
	else
	{
		G_SpawnString("angle", NOVALUE, &value);
		if (Q_stricmp(value, NOVALUE) != 0)
		{
			angles[YAW] = atof( value );
		}
		else
		{
			angles[YAW] = 0.0;
		}
		angles[YAW] = fmod(angles[YAW] + level.mRotationAdjust, 360.0f);
		Com_sprintf(temp, sizeof( temp ), "%0.0f", angles[YAW]);
		AddSpawnField("angle", temp);
	}

	// RJR experimental code for handling "direction" field of breakable brushes
	// though direction is rarely ever used.
	G_SpawnString("direction", NOVALUE, &value);
	if (Q_stricmp(value, NOVALUE) != 0)
	{
		if ( sscanf( value, "%f %f %f", &angles[0], &angles[1], &angles[2] ) != 3 ) {
			G_Printf( "HandleEntityAdjustment: failed sscanf on 'direction' (%s)\n", value );
			VectorClear( angles );
		}
	}
	else
	{
		angles[0] = angles[1] = angles[2] = 0.0;
	}
	angles[YAW] = fmod(angles[YAW] + level.mRotationAdjust, 360.0f);
	Com_sprintf(temp, sizeof( temp ), "%0.0f %0.0f %0.0f", angles[0], angles[1], angles[2]);
	AddSpawnField("direction", temp);


	AddSpawnField("BSPInstanceID", level.mTargetAdjust);

	G_SpawnString("targetname", NOVALUE, &value);
	if (Q_stricmp(value, NOVALUE) != 0)
	{
		Com_sprintf(temp, sizeof( temp ), "%s%s", level.mTargetAdjust, value);
		AddSpawnField("targetname", temp);
	}

	G_SpawnString("target", NOVALUE, &value);
	if (Q_stricmp(value, NOVALUE) != 0)
	{
		Com_sprintf(temp, sizeof( temp ), "%s%s", level.mTargetAdjust, value);
		AddSpawnField("target", temp);
	}

	G_SpawnString("killtarget", NOVALUE, &value);
	if (Q_stricmp(value, NOVALUE) != 0)
	{
		Com_sprintf(temp, sizeof( temp ), "%s%s", level.mTargetAdjust, value);
		AddSpawnField("killtarget", temp);
	}

//.........这里部分代码省略.........

	double n2;\
    double XS;\
    double YS;

#define PJ_LIB__
# include	"projects.h"
PROJ_HEAD(gstmerc, "Gauss-Schreiber Transverse Mercator (aka Gauss-Laborde Reunion)")
	"\n\tCyl, Sph&Ell\n\tlat_0= lon_0= k_0=";
FORWARD(s_forward); /* spheroid */
	double L, Ls, sinLs1, Ls1;
	L= P->n1*lp.lam;
    Ls= P->c+P->n1*log(pj_tsfn(-1.0*lp.phi,-1.0*sin(lp.phi),P->e));
    sinLs1= sin(L)/cosh(Ls);
    Ls1= log(pj_tsfn(-1.0*asin(sinLs1),0.0,0.0));
    xy.x= (P->XS + P->n2*Ls1)*P->ra;
    xy.y= (P->YS + P->n2*atan(sinh(Ls)/cos(L)))*P->ra;
    /*fprintf(stderr,"fwd:\nL      =%16.13f\nLs     =%16.13f\nLs1    =%16.13f\nLP(%16.13f,%16.13f)=XY(%16.4f,%16.4f)\n",L,Ls,Ls1,lp.lam+P->lam0,lp.phi,(xy.x*P->a + P->x0)*P->to_meter,(xy.y*P->a + P->y0)*P->to_meter);*/
	return (xy);
}
INVERSE(s_inverse); /* spheroid */
	double L, LC, sinC;
	L= atan(sinh((xy.x*P->a - P->XS)/P->n2)/cos((xy.y*P->a - P->YS)/P->n2));
    sinC= sin((xy.y*P->a - P->YS)/P->n2)/cosh((xy.x*P->a - P->XS)/P->n2);
    LC= log(pj_tsfn(-1.0*asin(sinC),0.0,0.0));
    lp.lam= L/P->n1;
    lp.phi= -1.0*pj_phi2(P->ctx, exp((LC-P->c)/P->n1),P->e);
    /*fprintf(stderr,"inv:\nL      =%16.13f\nsinC   =%16.13f\nLC     =%16.13f\nXY(%16.4f,%16.4f)=LP(%16.13f,%16.13f)\n",L,sinC,LC,((xy.x/P->ra)+P->x0)/P->to_meter,((xy.y/P->ra)+P->y0)/P->to_meter,lp.lam+P->lam0,lp.phi);*/
	return (lp);
}
FREEUP; if (P) pj_dalloc(P); }
ENTRY0(gstmerc)

char Ziehe_KreisNO90(short xr1,short zp1,short xr2,short zp2,short xrm,short zpm,char g03g02)
{
  double sx=(double)zp1,sy=(double)xr1,
	 ex=(double)zp2,ey=(double)xr2,
	 mx=(double)zpm,my=(double)xrm;
  double dx,dy,w,cosw,sinw;
  double rad1,rad2,stw,endw,step=(double)schritt_g02;
  char quadr1,quadr2,endprg=0;
  short x,y;

  dx=sx-mx;
  dy=sy-my;
  rad1=sqrt(dx*dx+dy*dy);
  cosw=dx/rad1;
  sinw=dy/rad1;
  stw=arc_wi(sinw,cosw,&quadr1);

  dx=ex-mx;
  dy=ey-my;
  rad2=sqrt(dx*dx+dy*dy);
  cosw=dx/rad2;
  sinw=dy/rad2;
  endw=arc_wi(sinw,cosw,&quadr2);

  if (quadr1!=quadr2) { /* Routine zZ noch nicht in Funktion */
			Meldung("M99 Quadranten-Fehler!");
			Status("Fehler!");
			getch();
			endprg=1;
		      }

  if (fabs(fabs(rad1)-fabs(rad2))>160)
    {
      Meldung("M99-0,8mm-Toleranz �berschritten");
      Status("Fehler");
      getch();
      endprg=1; /* Nur zZ */
    }

  if (!endprg)
    {
      if (g03g02)
	{ /* G02 */
	   if (endw<stw) stw-=360.0;
	   for(w=stw;w<endw && !endprg;w+=step)
	      {
		x=(short)(mx+rad1*cos(w*deg));
		y=(short)(my+rad1*sin(w*deg));
		Setzte_Stahl(y,x,tool);
		Koordinaten(y,x);
		xrad=y;
		zpos=x;
		Wait(prgconfig.wait_g01);
		endprg=Check_FlagTasten();
	      }
	} else { /* g03 */
		 if (endw>stw) endw-=360.0;
		 for(w=stw;w>endw && !endprg;w-=step)
		  {
		    x=(short)(mx+rad1*cos(w*deg));
		    y=(short)(my+rad1*sin(w*deg));
		    Setzte_Stahl(y,x,tool);
		    Koordinaten(y,x);
		    xrad=y;
		    zpos=x;
		    Wait(prgconfig.wait_g01);
		    endprg=Check_FlagTasten();
		  }
	       } /* von g03 */
    } /* von if !endprg */

    if (!endprg)
       {
	 Setzte_Stahl(xr2,zp2,tool);
	 xrad=xr2;
	 zpos=zp2;
	 Koordinaten(xr2,zp2);
       }

  return endprg;
}

int main(void) {
//	int i, j;
	int sine_table[NTABLE], cosine_table[NTABLE];
	int mseq_table[CYCLES * NTABLE * MSEQ_LENGTH];
	float size;
	float s;
//	float aI, aQ;
//	float uI[M];                // Store last M samples of yI
//	float uQ[M];                // Store last M samples of yQ
//	float vI, vQ;
//	float xI, xQ;
//	float w = 0;
//	float wq = 0;
	int spb = MSEQ_LENGTH * CYCLES * NTABLE;
	int last_samples[spb * MSG_LEN];
	int i = 0;

	USART_initialize(3000000, 1); 	// set baud rate to 3,000,000

	DAC_initialize(FS); 			// set sample rate
	DAC_ready_flag = 0; 		// clear the DAC ready flag
	DAC_start();				// start the DAC running
	ADC_initialize(FS);
	ADC_start();
	ADC_ready_flag = 0; 		// clear the ADC ready flag
	int mseq[MSEQ_LENGTH] = { 1, 0, 1, 0, 1, 0, 1 };
	int msg_table[MSG_LEN * CYCLES * NTABLE];
	int recieved[MSG_LEN];
	//char message[MSG_LEN] = "testing";
	int msg[MSG_LEN] = { 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0 };
	size = (float) (NTABLE / TABCYCLES);
	for (i = 0; i < NTABLE; i++) {
		sine_table[i] = (int) floor(32767 * sin(2.0 * PI * i / (float) size));
		cosine_table[i] = (int) floor(32767 * cos(2.0 * PI * i / (float) size));
	}
	for (i = 0; i < spb; i++) {
		mseq_table[i] =
				(mseq[i / (NTABLE * CYCLES)]) ?
						-sine_table[i % NTABLE] : cosine_table[i % NTABLE];
	}
	for (i = 0; i < CYCLES * NTABLE * MSG_LEN; i++) {
		msg_table[i] =
				(msg[i / (CYCLES * NTABLE)]) ?
						-sine_table[i % NTABLE] : sine_table[i % NTABLE];
	}
	for (i = 0; i < MSG_LEN; i++) {
		recieved[i] = 0;
	}

	i = 0;
	int m = 0;
	int strength = 0;
	int maxstrength = 0;
	for (i = 0; i < spb; i++) {
		last_samples[i] = 0;
	}
	int cycles = 0;
	int seconds = 0;
	int state = 0; // 0: listening for mseq
				   // 4: initializing collection of packet data
				   // 5: collecting packet data
				   // 6: decoding & printing data
	while (1) {

		if (ADC_ready_flag != 0 && DAC_ready_flag != 0) {
			ADC_ready_flag = 0;
			DAC_ready_flag = 0;

			if (state == 0) {
				DAC1_value = mseq_table[m];
				s = ADC1_value + 15000;                   // PIN 35A - Grab ADC value
				last_samples[m] = s;

				strength = 0;
				for (i = 0; i < spb; i++) { // Slow?
					strength +=
							(last_samples[(m + 1 + i) % spb] * mseq_table[i])
									/ 32767;
				}
				strength /= 100;
				m++;
				m = m % spb;
				cycles++;
				cycles = cycles % FS;
				if (strength > 7500)
					printf("Convolution: %i \n", strength);
				if (strength > 7500) {
					printf("CAUGHT MSEQ\n");
					state = 4;
					m = 0;
				}
				if (strength > maxstrength) {
					maxstrength = strength;
					printf("Max Strength: ");
					printf("%i", maxstrength);
					printf("\n");
				}
				if (cycles == 0) {
					seconds++;
					printf("[O3] Seconds elapsed: ");
//.........这里部分代码省略.........

StkFloat Whistle :: tick( unsigned int )
{
  StkFloat soundMix, tempFreq;
  StkFloat envOut = 0, temp, temp1, temp2, tempX, tempY;
  double phi, cosphi, sinphi;
  double gain = 0.5, mod = 0.0;

  if ( --subSampCount_ <= 0 )	{
    tempVectorP_ = pea_.getPosition();
    subSampCount_ = subSample_;
    temp = bumper_.isInside( tempVectorP_ );
#ifdef WHISTLE_ANIMATION
    frameCount += 1;
    if ( frameCount >= (1470 / subSample_) ) {
      frameCount = 0;
      printf("%f %f %f\n",tempVectorP_->getX(),tempVectorP_->getY(),envOut);
      fflush(stdout);
    }
#endif
    envOut = envelope_.tick();

    if (temp < (BUMP_RADIUS + PEA_RADIUS)) {
      tempX = envOut * tickSize_ * 2000 * noise_.tick();
      tempY = -envOut * tickSize_ * 1000 * (1.0 + noise_.tick());
      pea_.addVelocity( tempX, tempY, 0 ); 
      pea_.tick( tickSize_ );
    }

    mod  = exp(-temp * 0.01);	  // exp. distance falloff of fipple/pea effect
    temp = onepole_.tick(mod);	// smooth it a little
    gain = (1.0 - (fippleGainMod_*0.5)) + (2.0 * fippleGainMod_ * temp);
    gain *= gain;	              // squared distance/gain
    //    tempFreq = 1.0				//  Normalized Base Freq
    //			+ (fippleFreqMod_ * 0.25) - (fippleFreqMod_ * temp) // fippleModulation 
    //			- (blowFreqMod_) + (blowFreqMod_ * envOut); // blowingModulation
    // short form of above
    tempFreq = 1.0 + fippleFreqMod_*(0.25-temp) + blowFreqMod_*(envOut-1.0);
    tempFreq *= baseFrequency_;

    sine_.setFrequency(tempFreq);

    tempVectorP_ = pea_.getPosition();
    temp = can_.isInside(tempVectorP_);
    temp  = -temp;       // We know (hope) it's inside, just how much??
    if (temp < (PEA_RADIUS * 1.25)) {
      pea_.getVelocity( &tempVector_ );  // This is the can/pea collision
      tempX = tempVectorP_->getX();     // calculation.  Could probably
      tempY = tempVectorP_->getY();     // simplify using tables, etc.
      phi = -atan2(tempY,tempX);

      cosphi = cos(phi);
      sinphi = sin(phi);
      temp1 = (cosphi*tempVector_.getX()) - (sinphi*tempVector_.getY());
      temp2 = (sinphi*tempVector_.getX()) + (cosphi*tempVector_.getY());
      temp1 = -temp1;
      tempX = (cosphi*temp1) + (sinphi*temp2);
      tempY = (-sinphi*temp1) + (cosphi*temp2);
      pea_.setVelocity(tempX, tempY, 0);
      pea_.tick(tickSize_);
      pea_.setVelocity( tempX*canLoss_, tempY*canLoss_, 0 );
      pea_.tick(tickSize_);
    }

    temp = tempVectorP_->getLength();	
    if (temp > 0.01) {
      tempX = tempVectorP_->getX();
      tempY = tempVectorP_->getY();
      phi = atan2( tempY, tempX );
      phi += 0.3 * temp / CAN_RADIUS;
      cosphi = cos(phi);
      sinphi = sin(phi);
      tempX = 3.0 * temp * cosphi;
      tempY = 3.0 * temp * sinphi;
    }
    else {
      tempX = 0.0;
      tempY = 0.0;
    }

    temp = (0.9 + 0.1*subSample_*noise_.tick()) * envOut * 0.6 * tickSize_;
    pea_.addVelocity( temp * tempX, (temp*tempY) - (GRAVITY*tickSize_), 0 );
    pea_.tick( tickSize_ );

    // bumper_.tick(0.0);
  }

  temp = envOut * envOut * gain / 2;
  soundMix = temp * ( sine_.tick() + ( noiseGain_*noise_.tick() ) );
  lastFrame_[0] = 0.20 * soundMix; // should probably do one-zero filter here

  return lastFrame_[0];
}

void lsat_bandctes(lsat_data * lsat, int i, char method,
		   double percent, int dark, double rayleigh)
{
    double pi_d2, sin_e, cos_v;

    /* TAUv  = at. transmittance surface-sensor */
    /* TAUz  = at. transmittance sun-surface    */
    /* Edown = diffuse sky spectral irradiance  */
    double TAUv, TAUz, Edown;

    pi_d2 = (double)(PI * lsat->dist_es * lsat->dist_es);
    sin_e = (double)(sin(D2R * lsat->sun_elev));
    cos_v = (double)(cos(D2R * (lsat->number < 4 ? 9.2 : 8.2)));

    /** Global irradiance on the sensor.
	 * Radiance to reflectance coefficient, only NO thermal bands.
	 * K1 and K2 variables are also utilized as thermal constants
     */
    if (lsat->band[i].thermal == 0) {
	switch (method) {
	case DOS2:
	    {
		TAUv = 1.;
		TAUz = (lsat->band[i].wavemax < 1.) ? sin_e : 1.;
		Edown = 0.;
		break;
	    }
	case DOS2b:
	    {
		TAUv = (lsat->band[i].wavemax < 1.) ? cos_v : 1.;
		TAUz = (lsat->band[i].wavemax < 1.) ? sin_e : 1.;
		Edown = 0.;
		break;
	    }
	case DOS3:
	    {
		double t;

		t = 2. / (lsat->band[i].wavemax + lsat->band[i].wavemin);
		t = 0.008569 * t * t * t * t * (1 + 0.0113 * t * t + 0.000013 * t * t * t * t);
		TAUv = exp(-t / cos_v);
		TAUz = exp(-t / sin_e);
		Edown = rayleigh;
		break;
	    }
	case DOS4:
	    {
		double Ro =
		    (lsat->band[i].lmax - lsat->band[i].lmin) * (dark -
								 lsat->
								 band[i].
								 qcalmin) /
		    (lsat->band[i].qcalmax - lsat->band[i].qcalmin) +
		    lsat->band[i].lmin;
		double Tv = 1.;
		double Tz = 1.;
		double Lp = 0.;

		do {
		    TAUz = Tz;
		    TAUv = Tv;
		    Lp = Ro - percent * TAUv * (lsat->band[i].esun * sin_e * TAUz + PI * Lp) / pi_d2;
		    Tz = 1. - (4. * pi_d2 * Lp) / (lsat->band[i].esun * sin_e);
		    if (Tz <= 0)
			G_fatal_error(_("The DOS4 method is not applicable here: approximation of atmospheric transmittance coefficients is unstable. Use another DOS method or use other sun_elevation parameter"));
		    Tv = exp(sin_e * log(Tz) / cos_v);
		} while (TAUv != Tv && TAUz != Tz);
		TAUz = (Tz < 1. ? Tz : 1.);
		TAUv = (Tv < 1. ? Tv : 1.);
		Edown = (Lp < 0. ? 0. : PI * Lp);
		break;
	    }
	default:		/* DOS1 and Without atmospheric-correction */
	    TAUv = 1.;
	    TAUz = 1.;
	    Edown = 0.;
	    break;
	}
	lsat->band[i].K2 = 0.;
	lsat->band[i].K1 = TAUv * (lsat->band[i].esun * sin_e * TAUz + Edown) / pi_d2;
	if (method > DOS)
	    G_verbose_message("... TAUv = %.5f, TAUz = %.5f, Edown = %.5f\n", TAUv, TAUz, Edown);
    }

    /** Digital number to radiance coefficients.
	 * Without atmospheric calibration for thermal bands.
     */
    lsat->band[i].gain = (lsat->band[i].lmax - lsat->band[i].lmin) / (lsat->band[i].qcalmax - lsat->band[i].qcalmin);

    if (method == UNCORRECTED || lsat->band[i].thermal) {
	/* L = G * (DN - Qmin) + Lmin
	 *  -> bias = Lmin - G * Qmin    
	 */
	lsat->band[i].bias = (lsat->band[i].lmin - lsat->band[i].gain * lsat->band[i].qcalmin);
    }
    else if (method > DOS) {
	/* L = Lsat - Lpath = G * DNsat + B - (G *  + B - p * rad_sun) 
	 *   = G * DNsat - G *  + p * rad_sun
	 *  -> bias = p * rad_sun - G 
	 */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
		x->in2_connected = 0;
		x->in3_connected = 0;
		x->sync = 0;
		x->playthrough = 0;
		x->frame_increment = 1.0;
		x->verbose = 0;
		x->table = (float *) getbytes(x->table_length * sizeof(float));
		x->Wanal = (float *) getbytes(MAX_Nw*sizeof(float));	
		x->Wsyn = (float *) getbytes(MAX_Nw*sizeof(float));	
		x->input = (float *) getbytes(MAX_Nw*sizeof(float));	
		x->Hwin = (float *) getbytes(MAX_Nw*sizeof(float));
		x->bindex = (float *) getbytes( (MAX_N+1) * sizeof(float) );
		x->buffer = (float *) getbytes(MAX_N*sizeof(float));
		x->channel = (float *) getbytes((MAX_N+2)*sizeof(float));
		x->output = (float *) getbytes(MAX_Nw*sizeof(float));
		x->bitshuffle = (int *) getbytes((MAX_N*2)*sizeof(int));
		x->trigland = (float *) getbytes((MAX_N*2)*sizeof(float));
		x->c_lastphase_in = (float *) getbytes((MAX_N2+1)*sizeof(float));
		x->c_lastphase_out = (float *) getbytes((MAX_N2+1)*sizeof(float));
		x->lastamp = (float *) getbytes((MAX_N+1) * sizeof(float));
		x->lastfreq = (float *) getbytes((MAX_N+1) * sizeof(float));
		x->local_frame = (float *) getbytes((MAX_N+2)*sizeof(float));
		x->loveboat = (float **) getbytes(x->framecount*sizeof(float *));
		

		/* here we stay with old reallocation approach and pray */
		for(i=0;i<x->framecount;i++){
			x->loveboat[i] = (float *) getbytes(((x->N)+2)*sizeof(float));
			if(x->loveboat[i] == NULL){
				error("memory error");
				return;
			}
			memset((char *)x->loveboat[i],0,(x->N+2)*sizeof(float));
		}
	} else if(initialized == 1) {
		//free and allocate
		oldsize = (oldN+2)*sizeof(float);
		for(i = 0; i < last_framecount; i++){
			freebytes(x->loveboat[i],oldsize) ;
		}
		oldsize = last_framecount*sizeof(float *);
		freebytes(x->loveboat,oldsize);
		x->loveboat = (float **) getbytes(x->framecount*sizeof(float *));  
		for(i=0;i<x->framecount;i++){
			x->loveboat[i] = (float *) getbytes((x->N+2)*sizeof(float));
			if(x->loveboat[i] == NULL){
				error("memory error");
				return;
			}
			memset((char *)x->loveboat[i],0,(x->N+2)*sizeof(float));
		}
	}  
	memset((char *)x->input,0,x->Nw * sizeof(float));
	memset((char *)x->output,0,x->Nw * sizeof(float));
	memset((char *)x->c_lastphase_in,0,(x->N2+1) * sizeof(float));
	memset((char *)x->lastamp,0,(x->N+1)*sizeof(float));
	memset((char *)x->lastfreq,0,(x->N+1)*sizeof(float));
	memset((char *)x->bindex,0,(x->N+1)*sizeof(float));
	memset((char *)x->buffer,0,x->N * sizeof(float));
	if(!x->vector_size){
		post("zero vector size - something is really screwed up here!");
		return;
	}
	for ( i = 0; i < x->table_length; i++ ) {
		x->table[i] = (float) x->N * cos((float)i * TWOPI / (float)x->table_length);
	}	
	x->c_fundamental =  (float) x->R/(float)x->N ;
	x->c_factor_in =  (float) x->R/((float)x->vector_size * TWOPI);
	
	if( x->hi_freq < x->c_fundamental ) {
		x->hi_freq = x->topfreq ;
	}
	x->hi_bin = 1;  
	x->curfreq = 0;
	while( x->curfreq < x->hi_freq ) {
		++(x->hi_bin);
		x->curfreq += x->c_fundamental ;
	}
	
	x->lo_bin = 0;  
	x->curfreq = 0;
	while( x->curfreq < x->lo_freq ) {
		++(x->lo_bin);
		x->curfreq += x->c_fundamental ;
	}
	
	if( x->hi_bin > x->N2)
		x->hi_bin = x->N2 ;
	if(x->lo_bin > x->hi_bin)
		x->lo_bin = x->hi_bin;
	
	x->i_vector_size = 1.0/x->vector_size;
	x->pitch_increment = x->P*x->table_length/x->R;
	
	makewindows( x->Hwin, x->Wanal, x->Wsyn, x->Nw, x->N, x->D);
	init_rdft( x->N, x->bitshuffle, x->trigland);

		  
	x->lock = 0;
}

static void assign_param(t_functype ftype,t_iparams *newparam,
			 real old[MAXFORCEPARAM],int comb,double reppow)
{
  int  i,j;
  real tmp;

  /* Set to zero */
  for(j=0; (j<MAXFORCEPARAM); j++) 
    {
      newparam->generic.buf[j]=0.0;
    }
  switch (ftype) {
  case F_G96ANGLES:
    /* Post processing of input data: store cosine iso angle itself */
    newparam->harmonic.rA =cos(old[0]*DEG2RAD);
    newparam->harmonic.krA=old[1];
    newparam->harmonic.rB =cos(old[2]*DEG2RAD);
    newparam->harmonic.krB=old[3];
    break;
  case F_G96BONDS:
    /* Post processing of input data: store square of length itself */
    newparam->harmonic.rA =sqr(old[0]);
    newparam->harmonic.krA=old[1];
    newparam->harmonic.rB =sqr(old[2]);
    newparam->harmonic.krB=old[3];
    break;
  case F_FENEBONDS:
    newparam->fene.bm=old[0];
    newparam->fene.kb=old[1];
    break;
  case F_RESTRBONDS:
    newparam->restraint.lowA = old[0];
    newparam->restraint.up1A = old[1];
    newparam->restraint.up2A = old[2];
    newparam->restraint.kA   = old[3];
    newparam->restraint.lowB = old[4];
    newparam->restraint.up1B = old[5];
    newparam->restraint.up2B = old[6];
    newparam->restraint.kB   = old[7];
    break;
  case F_TABBONDS:
  case F_TABBONDSNC:
  case F_TABANGLES:
  case F_TABDIHS:
    newparam->tab.table = round_check(old[0],0,ftype,"table index");
    newparam->tab.kA    = old[1];
    newparam->tab.kB    = old[3];
    break;
  case F_CROSS_BOND_BONDS:
    newparam->cross_bb.r1e=old[0];
    newparam->cross_bb.r2e=old[1];
    newparam->cross_bb.krr=old[2];
    break;
  case F_CROSS_BOND_ANGLES:
    newparam->cross_ba.r1e=old[0];
    newparam->cross_ba.r2e=old[1];
    newparam->cross_ba.r3e=old[2];
    newparam->cross_ba.krt=old[3];
    break;
  case F_UREY_BRADLEY:
    newparam->u_b.thetaA=old[0];
    newparam->u_b.kthetaA=old[1];
    newparam->u_b.r13A=old[2];
    newparam->u_b.kUBA=old[3];
    newparam->u_b.thetaB=old[4];
    newparam->u_b.kthetaB=old[5];
    newparam->u_b.r13B=old[6];
    newparam->u_b.kUBB=old[7];
    break;
  case F_QUARTIC_ANGLES:
    newparam->qangle.theta=old[0];
    for(i=0; i<5; i++)
      newparam->qangle.c[i]=old[i+1];
    break;
  case F_LINEAR_ANGLES:
    newparam->linangle.aA    = old[0];
    newparam->linangle.klinA = old[1];
    newparam->linangle.aB    = old[2];
    newparam->linangle.klinB = old[3];
    break;
  case F_BONDS:
  case F_ANGLES:
  case F_HARMONIC:
  case F_IDIHS:
    newparam->harmonic.rA =old[0];
    newparam->harmonic.krA=old[1];
    newparam->harmonic.rB =old[2];
    newparam->harmonic.krB=old[3];
    break;
  case F_MORSE:
    newparam->morse.b0A    =old[0];
    newparam->morse.cbA    =old[1];
    newparam->morse.betaA  =old[2];
    newparam->morse.b0B    =old[3];
    newparam->morse.cbB    =old[4];
    newparam->morse.betaB  =old[5];
    break;
  case F_CUBICBONDS:
    newparam->cubic.b0    =old[0];
    newparam->cubic.kb    =old[1];
//.........这里部分代码省略.........

int  HBblocks(char * fname)
{ FILE * f=fopen(fname,"a");
  double tb,sb,cb,alpha,sa,ca,ta,samb,camb,dMb,MbHl,MbSM,MbH,MbH3,Q;
  if(!f) return 1;
  Q=findValW("Q");
  if(slhaDecayExists(pNum("h")) <0)  slhaDecayPrint("h", f);
  if(slhaDecayExists(pNum("H")) <0)  slhaDecayPrint("H", f);
  if(slhaDecayExists(pNum("H3"))<0)  slhaDecayPrint("H3",f);
  if(slhaDecayExists(pNum("t")) <0)  slhaDecayPrint("t", f);
  if(slhaDecayExists(pNum("H+"))<0)  slhaDecayPrint("H+",f);

  tb=findValW("tB");  
    sb=tb/sqrt(1+tb*tb);
    cb=1/sqrt(1+tb*tb);
  alpha=findValW("alpha");
    sa=sin(alpha);
    ca=cos(alpha);
    ta=sa/ca;
    samb=sa*cb-ca*sb;
    camb=ca*cb+sa*sb;
  dMb=findValW("dMb");


  MbSM=findValW("Mb");
  MbH= MbSM/(1+dMb)*(1+dMb*ta/tb);  
  MbH3=MbSM/(1+dMb)*(1-dMb/tb/tb);
  MbHl=MbSM/(1+dMb)*(1-dMb/ta/tb);
 
  fprintf(f,"Block HiggsBoundsInputHiggsCouplingsBosons\n");
  fprintf(f,"# Effective coupling normalised to SM one and squared\n");
  fprintf(f,"# For (*) normalized on Sin(2*W)\n"); 
  fprintf(f," %12.4E  3    25    24    24 # higgs-W-W \n",        SQR(samb)  );
  fprintf(f," %12.4E  3    25    23    23 # higgs-Z-Z \n",        SQR(samb)  );
  fprintf(f," %12.4E  3    25    25    23 # higgs-higgs-Z \n",    0.   );

  { assignVal("Q",pMass("h"));
    calcMainFunc();
    fprintf(f," %12.4E  3    25    21    21 # higgs-gluon-gluon\n",  SQR(findValW("LGGh")/findValW("LGGSM")) );           
    fprintf(f," %12.4E  3    25    22    22 # higgs-gamma-gamma\n",  SQR(findValW("LAAh")/findValW("LAASM")) );
  }                          
 
  fprintf(f," %12.4E  3    35    24    24 # higgs-W-W \n",        SQR(camb)  );
  fprintf(f," %12.4E  3    35    23    23 # higgs-Z-Z \n",        SQR(camb)  );
  fprintf(f," %12.4E  3    35    25    23 # higgs-higgs-Z \n",    0.  );
  fprintf(f," %12.4E  3    35    35    23 # higgs-higgs-Z \n",    0.  );
  
  
  { assignVal("Q",pMass("H"));
    calcMainFunc();
    fprintf(f," %12.4E  3    35    21    21 # higgs-gluon-gluon\n",SQR(findValW("LGGH")/findValW("LGGSM"))  );   
    fprintf(f," %12.4E  3    35    22    22 # higgs-gamma-gamma\n",SQR(findValW("LAAH")/findValW("LAASM"))  );
  }
  

  fprintf(f," %12.4E  3    36    24    24 # higgs-W-W \n",        0.  );
  fprintf(f," %12.4E  3    36    23    23 # higgs-Z-Z \n",        0.  );
  
  { assignVal("Q",pMass("H3"));
    calcMainFunc();
    fprintf(f," %12.4E  3    36    21    21 # higgs-gluon-gluon\n",SQR(findValW("LGGH3")/2/findValW("LGGSM")) );
    fprintf(f," %12.4E  3    36    22    22 # higgs-gamma-gamma\n",SQR(findValW("LAAH3")/2/findValW("LAASM")) );             
  }
  
  fprintf(f," %12.4E  3    36    25    23 #*higgs-higgs-Z \n",    SQR(camb)  );
  fprintf(f," %12.4E  3    36    35    23 #*higgs-higgs-Z \n",    SQR(samb)  );
  fprintf(f," %12.4E  3    36    36    23 #* higgs-higgs-Z \n",   0.  );

  fprintf(f,"Block HiggsBoundsInputHiggsCouplingsFermions\n");
  fprintf(f,"# Effective coupling normalised to SM one and squared\n");
  fprintf(f," %12.4E   %12.4E   3    25     5    5 # higgs-b-b \n"    ,SQR((sa/cb)*(MbHl/MbSM)),0.);
  fprintf(f," %12.4E   %12.4E   3    25     6    6 # higgs-top-top \n",SQR(ca/sb)              ,0.);
  fprintf(f," %12.4E   %12.4E   3    25    15   15 # higgs-tau-tau \n",SQR(sa/cb)              ,0.);

  fprintf(f," %12.4E   %12.4E   3    35     5    5 # higgs-b-b \n"    ,SQR((ca/cb)*(MbH/MbSM))  ,0.);
  fprintf(f," %12.4E   %12.4E   3    35     6    6 # higgs-top-top \n",SQR(sa/sb)              ,0.);  
  fprintf(f," %12.4E   %12.4E   3    35    15   15 # higgs-tau-tau \n",SQR(ca/cb)  ,0.);

  fprintf(f," %12.4E   %12.4E   3    36     5    5 # higgs-b-b \n"    ,0.,SQR(tb*(MbH3/MbSM)));
  fprintf(f," %12.4E   %12.4E   3    36     6    6 # higgs-top-top \n",0.,SQR(1/tb)          );
  fprintf(f," %12.4E   %12.4E   3    36    15   15 # higgs-tau-tau \n",0.,SQR(tb)            );
  
  assignValW("Q",Q);
  calcMainFunc();    
  fclose(f);
   
  return 0;
}

/*
  Program to test polygon functionality
*/
int main(int argc, char *argv[]) {
  Image *src;
  const int rows = 300;
  const int cols = 400;
  Polygon *p;
  Color Red;
  Color Orange;
  Color White;
  Color Blue;
  Point pt[100];
  int i;
  
  srand(42);

  color_set(&Red, 0.9, 0.2, 0.1 );
  color_set(&Orange, 0.95, 0.7, 0.3 );
  color_set(&White, 1.0, 1.0, 1.0 );
  color_set(&Blue, 0.2, 0.1, 0.95 );

  src = image_create(rows, cols);

  // make a simple square to test proper areas and locations
  // the square ought to be 20x20, include pixel (30,30) and exclude pixel (50, 50)
  point_set2D(&(pt[0]), 30, 30);
  point_set2D(&(pt[1]), 50, 30);
  point_set2D(&(pt[2]), 50, 50);
  point_set2D(&(pt[3]), 30, 50);

  p = polygon_createp(4, pt);

  printf("drawing a square\n");
  polygon_drawFill(p, src, Blue);

  // something more interesting
  for(i=0;i<50;i++) {
    float dr = rand() % 20;
    point_set2D(&(pt[i]), 
    200 + cos((float)i * M_PI * 2.0 / 50.0)*(70 + dr),
    150 + sin((float)i * M_PI * 2.0 / 50.0)*(70 + dr));
  }
  polygon_set(p, 50, pt);

  printf("drawing first big polygon\n");
  polygon_drawFill(p, src, Red);

  for(i=0;i<50;i++) {
    float dr = rand() % 15;
    point_set2D(&(pt[i]), 
    200 + cos((float)i * M_PI * 2.0 / 50.0)*(50 + dr),
    150 + sin((float)i * M_PI * 2.0 / 50.0)*(50 + dr));
  }
  polygon_set(p, 50, pt);

  printf("drawing second big polygon\n");
  polygon_drawFill(p, src, Orange);

  for(i=0;i<50;i++) {
    float dr = rand() % 10;
    point_set2D(&(pt[i]), 
    200 + cos((float)i * M_PI * 2.0 / 50.0)*(30 + dr),
    150 + sin((float)i * M_PI * 2.0 / 50.0)*(30 + dr));
  }
  polygon_set(p, 50, pt);

  printf("drawing third big polygon\n");
  polygon_drawFill(p, src, White);

  printf("writing output\n");
  image_write(src, "test4a.ppm");

  image_free(src);

  return(0);
}

static void *SpVali_RayPath_Thread(void *tid_p)
{
	size_t tid = *((size_t *)tid_p);
	Zone *zp, *zone;
	SpPhys *pp, *zone_pp;
	size_t i, plane, izone;
	GeRay ray;
	double r, theta, phi, t;
	gsl_rng *rng = gsl_rng_alloc(gsl_rng_ranlux);

	gsl_rng_set(rng, (unsigned long)time(NULL));

	for(izone = 0, zone = Zone_GetMinLeaf(v_raypath_model.grid); zone; izone++, zone = Zone_AscendTree(zone)) {
		zone_pp = zone->data;

		/* Skip empty and non-leaf zones */
		if(!zone_pp->non_empty_leaf)
			continue;

		if((izone % Sp_MPISIZE) != Sp_MPIRANK)
			continue;
			
		if(((izone / Sp_MPISIZE) % Sp_NTHREAD) != tid)
			continue;

		#define NRAY 5000

		/* This samples a random number uniformly in the
		 * interval [0, 1) */
		#define RAND()\
			gsl_rng_uniform(rng)
		/* This samples a random number uniformly in the
		 * interval (0, 1) */
		#define PRAND()\
			gsl_rng_uniform_pos(rng)

		for(i = 0; i < NRAY; i++) {
			if(1) {
				ray = GeRay_Rand(rng, &zone->voxel);
			}
			else {
				Mem_BZERO(&ray);

				/* Init random ray origin */
				ray.e = zone->voxel.cen;
				r = GeRay_E(ray, 0);
				theta = PRAND() * PI;
				phi = RAND() * TWOPI;

				/* Convert to rectangular coordinates */
				GeRay_E(ray, 0) = r * sin(theta) * cos(phi);
				GeRay_E(ray, 1) = r * sin(theta) * sin(phi);
				GeRay_E(ray, 2) = r * cos(theta);

				/* Set random ray direction: first obtain direction in spherical
				 * coordinates then convert to rectangular coordinates */

				/* theta must be sampled with PRAND() or there may be a
				 * bias towards the theta=0 direction */
				theta = asin(2.0 * RAND() - 1.0) + 0.5 * PI;

				/* phi must be sampled with RAND() to avoid insufficient
				 * samples in the phi=0 direction */
				phi = PRAND() * TWOPI;

				/* Convert to rectangular coordinates */
				GeRay_D(ray, 0) = sin(theta) * cos(phi);
				GeRay_D(ray, 1) = sin(theta) * sin(phi);
				GeRay_D(ray, 2) = cos(theta);
			}

			zp = zone;

			/* Calculate radiative transfer along this direction */
			/* Propagate the ray through the cloud until we've
			 * reached the edge */
			while(zp) {
				/* Calculate path to next boundary */
				GeRay_TraverseVoxel(&ray, &zp->voxel, &t, &plane);

				/* Pointer to physical parameters associated with this zone */
				pp = zp->data;

				/* Do radiative transfer only if there are molecules in
				 * this zone */
				if(pp->non_empty_leaf) {
					zone_pp->ds += t;
				}

				/* Propagate ray to next position */
				ray = GeRay_Inc(&ray, t);

				/* Get next zone to traverse to */
				zp = Zone_GetNext(zp, plane, &ray);
			}
		}
		zone_pp->ds /= (double)NRAY;
	}

	return NULL;
//.........这里部分代码省略.........

// Generates a normal random variable with 
// zero mean and stddev s
double gaussian_rand(double s)
{
  double U1 = uniform_rand(1);
  double U2 = uniform_rand(1);
  return s*sqrt(-2*log(U1))*cos(2*PI*U2);
}

void Mist::createRing(float x, float y, float radius, float width) {
  offsets.push_back(vertices.size());

  float resolution = 60.0;
  float angle = TWO_PI/resolution;
  float outer_radius = radius + width;

  for(int i = 0; i < resolution; ++i) {
    float c0 = cos( (i + 0) * angle);
    float s0 = sin( (i + 0) * angle);
    float c1 = cos( (i + 1) * angle);
    float s1 = sin( (i + 1) * angle);
    
    // positions
    vec3 pa(c0 * radius, s0 * radius, 0.0f);
    vec3 pb(c1 * radius, s1 * radius, 0.0f);
    vec3 pc(c1 * outer_radius, s1 * outer_radius, 0.0f);
    vec3 pd(c0 * outer_radius, s0 * outer_radius, 0.0f);

    // texcoords
    float u0 = float(i+0)/resolution;
    float u1 = float(i+1)/resolution;
    vec3 ta(u0, 0.0f, 1.0f);
    vec3 tb(u1, 0.0f, 1.0f);
    vec3 tc(u1, 1.0f, 1.0f);
    vec3 td(u0, 1.0f, 1.0f);

    // calculate distances from the corners to the centers
    vec3 intersection;
    if(!intersect(pa, pc, pb, pd, intersection)) {
      printf("The vertices of the dist do not intersect. Error.\n");
      ::exit(EXIT_FAILURE);
    }

    float d0 = length(pa - intersection);
    float d1 = length(pb - intersection);
    float d2 = length(pc - intersection);
    float d3 = length(pd - intersection);

    ta = ta * ((d0 + d2)/d2);
    tb = tb * ((d1 + d3)/d3);
    tc = tc * ((d2 + d0)/d0);
    td = td * ((d3 + d1)/d1);

    // store the vertices
    VertexPT3 a(pa,ta);
    VertexPT3 b(pb,tb);
    VertexPT3 c(pc,tc);
    VertexPT3 d(pd,td);

    vertices.push_back(a);
    vertices.push_back(b);
    vertices.push_back(c);

    vertices.push_back(a);
    vertices.push_back(c);
    vertices.push_back(d);
  }

  cou