自学教程：C++ tan函数代码示例

/** * /brief Horizontal move of point position * This function uses an algorithm for an oblate spheroid earth model. * The algorithm is described here:  * http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf */int nmea_move_horz_ellipsoid(    const nmeaPOS *start_pos,   /**< Start position in radians */    nmeaPOS *end_pos,           /**< (O) Result position in radians */    double azimuth,             /**< Azimuth in radians */    double distance,            /**< Distance (km) */    double *end_azimuth         /**< (O) Azimuth at end position in radians */    ){    /* Variables */    double f, a, b, sqr_a, sqr_b;    double phi1, tan_U1, sin_U1, cos_U1, s, alpha1, sin_alpha1, cos_alpha1;    double sigma1, sin_alpha, sqr_cos_alpha, sqr_u, A, B;    double sigma_initial, sigma, sigma_prev, sin_sigma, cos_sigma, cos_2_sigmam, sqr_cos_2_sigmam, delta_sigma;    int remaining_steps;    double tmp1, phi2, lambda, C, L;        /* Check input */    NMEA_ASSERT(start_pos != 0);    NMEA_ASSERT(end_pos != 0);        if (fabs(distance) < 1e-12)    { /* No move */        *end_pos = *start_pos;        if ( end_azimuth != 0 ) *end_azimuth = azimuth;        return ! (NMEA_POSIX(isnan)(end_pos->lat) || NMEA_POSIX(isnan)(end_pos->lon));    } /* No move */    /* Earth geometry */    f = NMEA_EARTH_FLATTENING;    a = NMEA_EARTH_SEMIMAJORAXIS_M;    b = (1 - f) * a;    sqr_a = a * a;    sqr_b = b * b;        /* Calculation */    phi1 = start_pos->lat;    tan_U1 = (1 - f) * tan(phi1);    cos_U1 = 1 / sqrt(1 + tan_U1 * tan_U1);    sin_U1 = tan_U1 * cos_U1;    s = distance;    alpha1 = azimuth;    sin_alpha1 = sin(alpha1);    cos_alpha1 = cos(alpha1);    sigma1 = atan2(tan_U1, cos_alpha1);    sin_alpha = cos_U1 * sin_alpha1;    sqr_cos_alpha = 1 - sin_alpha * sin_alpha;    sqr_u = sqr_cos_alpha * (sqr_a - sqr_b) / sqr_b;     A = 1 + sqr_u / 16384 * (4096 + sqr_u * (-768 + sqr_u * (320 - 175 * sqr_u)));    B = sqr_u / 1024 * (256 + sqr_u * (-128 + sqr_u * (74 - 47 * sqr_u)));        /* Initialize iteration */    sigma_initial = s / (b * A);    sigma = sigma_initial;    sin_sigma = sin(sigma);    cos_sigma = cos(sigma);    cos_2_sigmam = cos(2 * sigma1 + sigma);    sqr_cos_2_sigmam = cos_2_sigmam * cos_2_sigmam;    delta_sigma = 0;    sigma_prev = 2 * NMEA_PI;    remaining_steps = 20;    while ((fabs(sigma - sigma_prev) > 1e-12) && (remaining_steps > 0))    { /* Iterate */        cos_2_sigmam = cos(2 * sigma1 + sigma);        sqr_cos_2_sigmam = cos_2_sigmam * cos_2_sigmam;        sin_sigma = sin(sigma);        cos_sigma = cos(sigma);        delta_sigma = B * sin_sigma * (              cos_2_sigmam + B / 4 * (              cos_sigma * (-1 + 2 * sqr_cos_2_sigmam) -              B / 6 * cos_2_sigmam * (-3 + 4 * sin_sigma * sin_sigma) * (-3 + 4 * sqr_cos_2_sigmam)             ));        sigma_prev = sigma;        sigma = sigma_initial + delta_sigma;        remaining_steps --;    } /* Iterate */        /* Calculate result */    tmp1 = (sin_U1 * sin_sigma - cos_U1 * cos_sigma * cos_alpha1);    phi2 = atan2(            sin_U1 * cos_sigma + cos_U1 * sin_sigma * cos_alpha1,            (1 - f) * sqrt(sin_alpha * sin_alpha + tmp1 * tmp1)            );    lambda = atan2(            sin_sigma * sin_alpha1,            cos_U1 * cos_sigma - sin_U1 * sin_sigma * cos_alpha1            );    C = f / 16 * sqr_cos_alpha * (4 + f * (4 - 3 * sqr_cos_alpha));    L = lambda -        (1 - C) * f * sin_alpha * (        sigma + C * sin_sigma *        (cos_2_sigmam + C * cos_sigma * (-1 + 2 * sqr_cos_2_sigmam))        );    //.........这里部分代码省略.........

/*===============CG_CalcViewValuesSets cg.refdef view values===============*/static int CG_CalcViewValues( void ) {	playerState_t   *ps;	static vec3_t oldOrigin = {0,0,0};	memset( &cg.refdef, 0, sizeof( cg.refdef ) );	// strings for in game rendering	// Q_strncpyz( cg.refdef.text[0], "Park Ranger", sizeof(cg.refdef.text[0]) );	// Q_strncpyz( cg.refdef.text[1], "19", sizeof(cg.refdef.text[1]) );	// calculate size of 3D view	CG_CalcVrect();	ps = &cg.predictedPlayerState;	if ( cg.cameraMode ) {		vec3_t origin, angles;		float fov = 90;		float x;		if ( trap_getCameraInfo( CAM_PRIMARY, cg.time, &origin, &angles, &fov ) ) {			VectorCopy( origin, cg.refdef.vieworg );			angles[ROLL] = 0;			angles[PITCH] = -angles[PITCH];     // (SA) compensate for reversed pitch (this makes the game match the editor, however I'm guessing the real fix is to be done there)			VectorCopy( angles, cg.refdefViewAngles );			AnglesToAxis( cg.refdefViewAngles, cg.refdef.viewaxis );			x = cg.refdef.width / tan( fov / 360 * M_PI );			cg.refdef.fov_y = atan2( cg.refdef.height, x );			cg.refdef.fov_y = cg.refdef.fov_y * 360 / M_PI;			cg.refdef.fov_x = fov;			// RF, had to disable, sometimes a loadgame to a camera in the same position			// can cause the game to not know where the camera is, therefore snapshots			// dont show the correct entities			//if(VectorCompare(origin, oldOrigin))			//	return 0;			VectorCopy( origin, oldOrigin );			trap_SendClientCommand( va( "setCameraOrigin %f %f %f", origin[0], origin[1], origin[2] ) );			return 0;		} else {			cg.cameraMode = qfalse;                 // camera off in cgame			trap_Cvar_Set( "cg_letterbox", "0" );			trap_SendClientCommand( "stopCamera" );    // camera off in game			trap_stopCamera( CAM_PRIMARY );           // camera off in client			CG_Fade( 0, 0, 0, 255, 0, 0 );                // go black			CG_Fade( 0, 0, 0, 0, cg.time + 200, 1500 );   // then fadeup		}	}	// intermission view	if ( ps->pm_type == PM_INTERMISSION ) {		VectorCopy( ps->origin, cg.refdef.vieworg );		VectorCopy( ps->viewangles, cg.refdefViewAngles );		AnglesToAxis( cg.refdefViewAngles, cg.refdef.viewaxis );		return CG_CalcFov();	}	cg.bobcycle = ( ps->bobCycle & 128 ) >> 7;	cg.bobfracsin = fabs( sin( ( ps->bobCycle & 127 ) / 127.0 * M_PI ) );	cg.xyspeed = sqrt( ps->velocity[0] * ps->velocity[0] +					   ps->velocity[1] * ps->velocity[1] );	VectorCopy( ps->origin, cg.refdef.vieworg );	VectorCopy( ps->viewangles, cg.refdefViewAngles );	// add error decay	if ( cg_errorDecay.value > 0 ) {		int t;		float f;		t = cg.time - cg.predictedErrorTime;		f = ( cg_errorDecay.value - t ) / cg_errorDecay.value;		if ( f > 0 && f < 1 ) {			VectorMA( cg.refdef.vieworg, f, cg.predictedError, cg.refdef.vieworg );		} else {			cg.predictedErrorTime = 0;		}	}	// Ridah, lock the viewangles if the game has told us to	if ( ps->viewlocked ) {		/*		if (ps->viewlocked == 4)		{			centity_t *tent;			tent = &cg_entities[ps->viewlocked_entNum];			VectorCopy (tent->currentState.apos.trBase, cg.refdefViewAngles);//.........这里部分代码省略.........

static int CG_CalcFov( void ) {	static float lastfov = 90;      // for transitions back from zoomed in modes	float x;	float phase;	float v;	int contents;	float fov_x, fov_y;	float zoomFov;	float f;	int inwater;	qboolean dead;	CG_Zoom();	if ( cg.predictedPlayerState.stats[STAT_HEALTH] <= 0 ) {		dead = qtrue;		cg.zoomedBinoc = qfalse;		cg.zoomTime = 0;		cg.zoomval = 0;	} else {		dead = qfalse;	}	if ( cg.predictedPlayerState.pm_type == PM_INTERMISSION ) {		// if in intermission, use a fixed value		fov_x = 90;	} else {		// user selectable		if ( cgs.dmflags & DF_FIXED_FOV ) {			// dmflag to prevent wide fov for all clients			fov_x = 90;		} else {			fov_x = cg_fov.value;			if ( fov_x < 1 ) {				fov_x = 1;			} else if ( fov_x > 160 ) {				fov_x = 160;			}		}		// account for zooms		if ( cg.zoomval ) {			zoomFov = cg.zoomval;   // (SA) use user scrolled amount			if ( zoomFov < 1 ) {				zoomFov = 1;			} else if ( zoomFov > 160 ) {				zoomFov = 160;			}		} else {			zoomFov = lastfov;		}		// do smooth transitions for the binocs		if ( cg.zoomedBinoc ) {        // binoc zooming in			f = ( cg.time - cg.zoomTime ) / (float)ZOOM_TIME;			if ( f > 1.0 ) {				fov_x = zoomFov;			} else {				fov_x = fov_x + f * ( zoomFov - fov_x );			}			lastfov = fov_x;		} else if ( cg.zoomval ) {    // zoomed by sniper/snooper			fov_x = cg.zoomval;			lastfov = fov_x;		} else {                    // binoc zooming out			f = ( cg.time - cg.zoomTime ) / (float)ZOOM_TIME;			if ( f > 1.0 ) {				fov_x = fov_x;			} else {				fov_x = zoomFov + f * ( fov_x - zoomFov );			}		}	}	// DHM - Nerve :: zoom in for Limbo or Spectator	if ( cgs.gametype == GT_WOLF ) {		if ( cg.snap->ps.pm_flags & PMF_FOLLOW && cg.snap->ps.weapon == WP_SNIPERRIFLE ) {			fov_x = cg_zoomDefaultSniper.value;		}	}	// dhm - end	if ( !dead && ( cg.weaponSelect == WP_SNOOPERSCOPE ) ) {		cg.refdef.rdflags |= RDF_SNOOPERVIEW;	} else {		cg.refdef.rdflags &= ~RDF_SNOOPERVIEW;	}	if ( cg.snap->ps.persistant[PERS_HWEAPON_USE] ) {		fov_x = 55;	}	x = cg.refdef.width / tan( fov_x / 360 * M_PI );	fov_y = atan2( cg.refdef.height, x );	fov_y = fov_y * 360 / M_PI;	// warp if underwater	contents = CG_PointContents( cg.refdef.vieworg, -1 );	if ( contents & ( CONTENTS_WATER | CONTENTS_SLIME | CONTENTS_LAVA ) ) {//.........这里部分代码省略.........

//@Jkent: Is there a way to get access to R_customheight/width and/or r_height/width? (Get the newer height-dependant FOV rather than width-dependant as 1.1 widescreens actually have smaller FOVs)static int CG_CalcFov( void ){  float     x;  float     phase;  float     v;  int       contents;  float     fov_x, fov_y;  float     zoomFov;  float     f;  int       inwater;  int       attribFov;  usercmd_t cmd;  int       cmdNum;  cmdNum = trap_GetCurrentCmdNumber( );  trap_GetUserCmd( cmdNum, &cmd );  if( cg.predictedPlayerState.pm_type == PM_INTERMISSION ||      ( cg.snap->ps.persistant[ PERS_TEAM ] == TEAM_SPECTATOR ) )  {    // if in intermission, use a fixed value    fov_x = 90;  }  else  {    // don't lock the fov globally - we need to be able to change it    attribFov = BG_FindFovForClass( cg.predictedPlayerState.stats[ STAT_PCLASS ] );    fov_x = attribFov;    if ( fov_x < 1 )      fov_x = 1;    else if ( fov_x > 160 )      fov_x = 160;    if( cg.spawnTime > ( cg.time - FOVWARPTIME ) &&        BG_ClassHasAbility( cg.predictedPlayerState.stats[ STAT_PCLASS ], SCA_FOVWARPS ) )    {      float temp, temp2;      temp = (float)( cg.time - cg.spawnTime ) / FOVWARPTIME;      temp2 = ( 170 - fov_x ) * temp;      //Com_Printf( "%f %f/n", temp*100, temp2*100 );      fov_x = 170 - temp2;    }    // account for zooms    zoomFov = BG_FindZoomFovForWeapon( cg.predictedPlayerState.weapon );    if ( zoomFov < 1 )      zoomFov = 1;    else if ( zoomFov > attribFov )      zoomFov = attribFov;    // only do all the zoom stuff if the client CAN zoom    // FIXME: zoom control is currently hard coded to BUTTON_ATTACK2    if( BG_WeaponCanZoom( cg.predictedPlayerState.weapon ) )    {      if ( cg.zoomed )      {        f = ( cg.time - cg.zoomTime ) / (float)ZOOM_TIME;        if ( f > 1.0 )          fov_x = zoomFov;        else          fov_x = fov_x + f * ( zoomFov - fov_x );        // BUTTON_ATTACK2 isn't held so unzoom next time        if( !( cmd.buttons & BUTTON_ATTACK2 ) )        {          cg.zoomed   = qfalse;          cg.zoomTime = cg.time;        }      }      else      {        f = ( cg.time - cg.zoomTime ) / (float)ZOOM_TIME;        if ( f > 1.0 )          fov_x = fov_x;        else          fov_x = zoomFov + f * ( fov_x - zoomFov );        // BUTTON_ATTACK2 is held so zoom next time        if( cmd.buttons & BUTTON_ATTACK2 )        {          cg.zoomed   = qtrue;          cg.zoomTime = cg.time;        }      }    }  }  x = cg.refdef.width / tan( fov_x / 360 * M_PI );  fov_y = atan2( cg.refdef.height, x );  fov_y = fov_y * 360 / M_PI;  // warp if underwater  contents = CG_PointContents( cg.refdef.vieworg, -1 );//.........这里部分代码省略.........

static intmsProjectShapeLine(projectionObj *in, projectionObj *out,                   shapeObj *shape, int line_index){  int i;  pointObj  lastPoint, thisPoint, wrkPoint, firstPoint;  lineObj *line = shape->line + line_index;  lineObj *line_out = line;  int valid_flag = 0; /* 1=true, -1=false, 0=unknown */  int numpoints_in = line->numpoints;  int line_alloc = numpoints_in;  int wrap_test;#ifdef USE_PROJ_FASTPATHS#define MAXEXTENT 20037508.34#define M_PIby360 .0087266462599716479#define MAXEXTENTby180 111319.4907777777777777777  if(in->wellknownprojection == wkp_lonlat && out->wellknownprojection == wkp_gmerc) {    for( i = line->numpoints-1; i >= 0; i-- ) {#define p_x line->point[i].x#define p_y line->point[i].y      p_x *= MAXEXTENTby180;      p_y = log(tan((90 + p_y) * M_PIby360)) * MS_RAD_TO_DEG;      p_y *= MAXEXTENTby180;      if (p_x > MAXEXTENT) p_x = MAXEXTENT;      if (p_x < -MAXEXTENT) p_x = -MAXEXTENT;      if (p_y > MAXEXTENT) p_y = MAXEXTENT;      if (p_y < -MAXEXTENT) p_y = -MAXEXTENT;#undef p_x#undef p_y    }    return MS_SUCCESS;  }#endif  wrap_test = out != NULL && out->proj != NULL && pj_is_latlong(out->proj)              && !pj_is_latlong(in->proj);  line->numpoints = 0;  if( numpoints_in > 0 )    firstPoint = line->point[0];  memset( &lastPoint, 0, sizeof(lastPoint) );  /* -------------------------------------------------------------------- */  /*      Loop over all input points in linestring.                       */  /* -------------------------------------------------------------------- */  for( i=0; i < numpoints_in; i++ ) {    int ms_err;    wrkPoint = thisPoint = line->point[i];    ms_err = msProjectPoint(in, out, &wrkPoint );    /* -------------------------------------------------------------------- */    /*      Apply wrap logic.                                               */    /* -------------------------------------------------------------------- */    if( wrap_test && i > 0 && ms_err != MS_FAILURE ) {      double dist;      pointObj pt1Geo;      if( line_out->numpoints > 0 )        pt1Geo = line_out->point[0];      else        pt1Geo = wrkPoint; /* this is a cop out */      dist = wrkPoint.x - pt1Geo.x;      if( fabs(dist) > 180.0          && msTestNeedWrap( thisPoint, firstPoint,                             pt1Geo, in, out ) ) {        if( dist > 0.0 )          wrkPoint.x -= 360.0;        else if( dist < 0.0 )          wrkPoint.x += 360.0;      }    }    /* -------------------------------------------------------------------- */    /*      Put result into output line with appropriate logic for          */    /*      failure breaking lines, etc.                                    */    /* -------------------------------------------------------------------- */    if( ms_err == MS_FAILURE ) {      /* We have started out invalid */      if( i == 0 ) {        valid_flag = -1;      }      /* valid data has ended, we need to work out the horizon */      else if( valid_flag == 1 ) {        pointObj startPoint, endPoint;        startPoint = lastPoint;        endPoint = thisPoint;        if( msProjectSegment( in, out, &startPoint, &endPoint )            == MS_SUCCESS ) {          line_out->point[line_out->numpoints++] = endPoint;        }//.........这里部分代码省略.........

double QgsDistanceArea::computeDistanceBearing(  const QgsPoint& p1, const QgsPoint& p2,  double* course1, double* course2 ){  if ( p1.x() == p2.x() && p1.y() == p2.y() )    return 0;  // ellipsoid  double a = mSemiMajor;  double b = mSemiMinor;  double f = 1 / mInvFlattening;  double p1_lat = DEG2RAD( p1.y() ), p1_lon = DEG2RAD( p1.x() );  double p2_lat = DEG2RAD( p2.y() ), p2_lon = DEG2RAD( p2.x() );  double L = p2_lon - p1_lon;  double U1 = atan(( 1 - f ) * tan( p1_lat ) );  double U2 = atan(( 1 - f ) * tan( p2_lat ) );  double sinU1 = sin( U1 ), cosU1 = cos( U1 );  double sinU2 = sin( U2 ), cosU2 = cos( U2 );  double lambda = L;  double lambdaP = 2 * M_PI;  double sinLambda = 0;  double cosLambda = 0;  double sinSigma = 0;  double cosSigma = 0;  double sigma = 0;  double alpha = 0;  double cosSqAlpha = 0;  double cos2SigmaM = 0;  double C = 0;  double tu1 = 0;  double tu2 = 0;  int iterLimit = 20;  while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )  {    sinLambda = sin( lambda );    cosLambda = cos( lambda );    tu1 = ( cosU2 * sinLambda );    tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );    sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );    cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;    sigma = atan2( sinSigma, cosSigma );    alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );    cosSqAlpha = cos( alpha ) * cos( alpha );    cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;    C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );    lambdaP = lambda;    lambda = L + ( 1 - C ) * f * sin( alpha ) *             ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );  }  if ( iterLimit == 0 )    return -1;  // formula failed to converge  double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );  double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );  double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );  double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -                                       B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );  double s = b * A * ( sigma - deltaSigma );  if ( course1 )  {    *course1 = atan2( tu1, tu2 );  }  if ( course2 )  {    // PI is added to return azimuth from P2 to P1    *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) + M_PI;  }  return s;}

	//function to find the body state vector at epoch	int body::locate_body(const double& epoch, double* state, const bool& need_deriv, missionoptions* options)	{		double DT, n, M, E, V[6];		switch (body_ephemeris_source)		{			case 1: //SPICE				double LT_dump;				spkez_c(spice_ID, epoch - (51544.5 * 86400.0), "J2000", "NONE", this->central_body_spice_ID, state, &LT_dump);				if (need_deriv)				{					double statepert[6];					spkez_c(spice_ID, epoch - (51544.5 * 86400.0) + 10.0, "J2000", "NONE", this->central_body_spice_ID, statepert, &LT_dump);					state[6] = (statepert[3] - state[3]) / (10.0);					state[7] = (statepert[4] - state[4]) / (10.0);					state[8] = (statepert[5] - state[5]) / (10.0);				}				break;			case 0: //static ephemeris					//TODO static ephemeris is not ready!					//note, always should give in Earth equatorial J2000 coordinates for internal processing					DT = ( epoch - this->reference_epoch );										if (this->SMA > 0.0)						n = sqrt(this->universe_mu / (this->SMA*this->SMA*this->SMA));					else						n = sqrt(this->universe_mu / (-this->SMA*this->SMA*this->SMA));										M = this->MA + n*DT;					M = fmod(M, 2 * EMTG::math::PI);					E = Kepler::KeplerLaguerreConway(this->ECC, M);					V[0] = this->SMA; 					V[1] = this->ECC;					V[2] = this->INC;					V[3] = this->RAAN;					V[4] = this->AOP;					true_anomaly = 2.0*atan(sqrt((1.0 + this->ECC) / (1.0 - this->ECC))*tan(E / 2.0));					V[5] = true_anomaly;					COE2inertial(V, this->universe_mu, state);					if (need_deriv)					{						double r = sqrt(state[0]*state[0] + state[1]*state[1] + state[2]*state[2]);						double r3 = r*r*r;						state[6] = -universe_mu/r3 * state[0];						state[7] = -universe_mu/r3 * state[1];						state[8] = -universe_mu/r3 * state[2];					}				break;			default:				cout << "Invalid ephemeris source " << body_ephemeris_source << " for object " << name << endl;				cout << "Program halted. Press enter to quit." << endl;#ifndef BACKGROUND_MODE				cin.ignore();#endif		}		return 0;	}

// Draws the FBO texture for 3DTV.void ApplicationOverlay::displayOverlayTexture3DTV(Camera& whichCamera, float aspectRatio, float fov) {    if (_alpha == 0.0f) {        return;    }        Application* application = Application::getInstance();        MyAvatar* myAvatar = application->getAvatar();    const glm::vec3& viewMatrixTranslation = application->getViewMatrixTranslation();        glActiveTexture(GL_TEXTURE0);        glEnable(GL_BLEND);    glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, GL_CONSTANT_ALPHA, GL_ONE);    _overlays.bindTexture();    glEnable(GL_DEPTH_TEST);    glDisable(GL_LIGHTING);    glEnable(GL_TEXTURE_2D);        glMatrixMode(GL_MODELVIEW);        glPushMatrix();    glLoadIdentity();    // Transform to world space    glm::quat rotation = whichCamera.getRotation();    glm::vec3 axis2 = glm::axis(rotation);    glRotatef(-glm::degrees(glm::angle(rotation)), axis2.x, axis2.y, axis2.z);    glTranslatef(viewMatrixTranslation.x, viewMatrixTranslation.y, viewMatrixTranslation.z);        // Translate to the front of the camera    glm::vec3 pos = whichCamera.getPosition();    glm::quat rot = myAvatar->getOrientation();    glm::vec3 axis = glm::axis(rot);        glTranslatef(pos.x, pos.y, pos.z);    glRotatef(glm::degrees(glm::angle(rot)), axis.x, axis.y, axis.z);        glColor4f(1.0f, 1.0f, 1.0f, _alpha);        //Render    const GLfloat distance = 1.0f;        const GLfloat halfQuadHeight = distance * tan(fov);    const GLfloat halfQuadWidth = halfQuadHeight * aspectRatio;    const GLfloat quadWidth = halfQuadWidth * 2.0f;    const GLfloat quadHeight = halfQuadHeight * 2.0f;        GLfloat x = -halfQuadWidth;    GLfloat y = -halfQuadHeight;    glDisable(GL_DEPTH_TEST);        glBegin(GL_QUADS);        glTexCoord2f(0.0f, 1.0f); glVertex3f(x, y + quadHeight, -distance);    glTexCoord2f(1.0f, 1.0f); glVertex3f(x + quadWidth, y + quadHeight, -distance);    glTexCoord2f(1.0f, 0.0f); glVertex3f(x + quadWidth, y, -distance);    glTexCoord2f(0.0f, 0.0f); glVertex3f(x, y, -distance);        glEnd();        if (_crosshairTexture == 0) {        _crosshairTexture = Application::getInstance()->getGLWidget()->bindTexture(QImage(Application::resourcesPath() + "images/sixense-reticle.png"));    }        //draw the mouse pointer    glBindTexture(GL_TEXTURE_2D, _crosshairTexture);        const float reticleSize = 40.0f / application->getGLWidget()->width() * quadWidth;    x -= reticleSize / 2.0f;    y += reticleSize / 2.0f;    const float mouseX = (application->getMouseX() / (float)application->getGLWidget()->width()) * quadWidth;    const float mouseY = (1.0 - (application->getMouseY() / (float)application->getGLWidget()->height())) * quadHeight;        glBegin(GL_QUADS);        glColor3f(RETICLE_COLOR[0], RETICLE_COLOR[1], RETICLE_COLOR[2]);        glTexCoord2d(0.0f, 0.0f); glVertex3f(x + mouseX, y + mouseY, -distance);    glTexCoord2d(1.0f, 0.0f); glVertex3f(x + mouseX + reticleSize, y + mouseY, -distance);    glTexCoord2d(1.0f, 1.0f); glVertex3f(x + mouseX + reticleSize, y + mouseY - reticleSize, -distance);    glTexCoord2d(0.0f, 1.0f); glVertex3f(x + mouseX, y + mouseY - reticleSize, -distance);        glEnd();        glEnable(GL_DEPTH_TEST);        glPopMatrix();        glDepthMask(GL_TRUE);    glBindTexture(GL_TEXTURE_2D, 0);    glDisable(GL_TEXTURE_2D);        glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, GL_CONSTANT_ALPHA, GL_ONE);    glEnable(GL_LIGHTING);        glColor4f(1.0f, 1.0f, 1.0f, 1.0f);}

treetan (tree t) {  if (is_double (t)) return as_tree (tan (as_double (t)));  return tree (TAN, t);}

float CalcFrustumScale(float fFovDeg) {	const float degToRad = 3.14159f * 2.0f / 360.0f;	float fFovRad = fFovDeg * degToRad;	return 1.0f / tan(fFovRad / 2.0f);}

void Raytracer::render(const char *filename, const char *depth_filename,                       Scene const &scene){    // Allocate the two images that will ultimately be saved.    Image colorImage(scene.resolution[0], scene.resolution[1]);    Image depthImage(scene.resolution[0], scene.resolution[1]);        // Create the zBuffer.    double *zBuffer = new double[scene.resolution[0] * scene.resolution[1]];    for(int i = 0; i < scene.resolution[0] * scene.resolution[1]; i++) {        zBuffer[i] = DBL_MAX;    }	// @@@@@@ YOUR CODE HERE	// calculate camera parameters for rays, refer to the slides for details	//!!! USEFUL NOTES: tan() takes rad rather than degree, use deg2rad() to transform	//!!! USEFUL NOTES: view plane can be anywhere, but it will be implemented differently,	//you can find references from the course slides 22_GlobalIllum.pdf        	Vector cameraPos = scene.camera.position;	Vector cameraCenter = scene.camera.center;		Vector cameraPosR = scene.camera.position;	Vector cameraCenterR = scene.camera.center;		// viewing direction vector get by taking center and subtracting camera position	Vector wVecOriginal = scene.camera.center; - cameraPos;    wVecOriginal.normalize();	// up vector is defined (u)    Vector uVec = scene.camera.up;    uVec.normalize();	// right vector is gotten by taking the cross product of w and v	Vector rVecOriginal = wVecOriginal.cross(uVec);	rVecOriginal.normalize();		double stereoDisplacement = scene.camera.stereoDist / 2.0;	int widthResolution = scene.resolution[0];	if (scene.camera.stereoDist > 0.0) {				printf("Start left picture./n");		cameraPos = scene.camera.position + (rVecOriginal * stereoDisplacement);		cameraPosR = scene.camera.position - (rVecOriginal * stereoDisplacement);				widthResolution = floor(scene.resolution[0] / 2);	} else if (scene.camera.stereoDist < 0.0) {				printf("Start left picture./n");		stereoDisplacement = - scene.camera.stereoDist / 2.0;		cameraPos = scene.camera.position - (rVecOriginal * stereoDisplacement);		cameraPosR = scene.camera.position + (rVecOriginal * stereoDisplacement);				widthResolution = floor(scene.resolution[0] / 2);	}	    Vector wVec = cameraCenter - cameraPos;    wVec.normalize();    Vector rVec = wVec.cross(uVec);    rVec.normalize();        // get top from tan(fovy)    double tangent = tan(deg2rad(scene.camera.fovy/2));	//double atangent = atan(deg2rad(scene.camera.fovy)/2);    // get length of top from centre of image plane    double top = scene.camera.zNear * tangent;    double right = top * scene.camera.aspect;	if (scene.camera.stereoDist != 0.0) {				right = right / 2;	}    double left = -right;    double bottom = -top;	    // calculate vector from camera to left top of image plane    Vector centerVec = cameraPos + (scene.camera.zNear * wVec);    Vector oVec = centerVec + (left * rVec) + (bottom * uVec);    double deltaU = (right - left) / scene.resolution[0];	if (scene.camera.stereoDist != 0.0) {				deltaU = deltaU * 2;	}    double deltaV = (top - bottom) / scene.resolution[1];    	        // Iterate over all the pixels in the image.    for(int y = 0; y < scene.resolution[1]; y++) {        for(int x = 0; x < widthResolution; x++) {            // Generate the appropriate ray for this pixel			Ray ray;			if (scene.objects.empty())			{				//no objects in the scene, then we render the default scene:				//in the default scene, we assume the view plane is at z = 640 with width and height both 640				ray = Ray(cameraPos, (Vector(-320, -320, 640) + Vector(x + 0.5, y + 0.5, 0) - cameraPos).normalized());			}			else			{				// set primary ray using the camera parameters				//!!! USEFUL NOTES: all world coordinate rays need to have a normalized direction								Vector changeU = (x + 0.5) * deltaU * rVec;                Vector changeY = (y + 0.5) * deltaV * uVec;                Vector pixelPos = oVec + changeU + changeY;//.........这里部分代码省略.........

/*=================UI_PlayerSetup_CalcFovassume refdef is valid=================*/static void UI_PlayerSetup_CalcFov( ref_params_t *fd ){	float x = fd->viewport[2] / tan( DEG2RAD( fd->fov_x ) * 0.5f );	float half_fov_y = atan( fd->viewport[3] / x );	fd->fov_y = RAD2DEG( half_fov_y ) * 2;}

//.........这里部分代码省略.........	uiPlayerSetup.topColor.generic.id = ID_TOPCOLOR;	uiPlayerSetup.topColor.generic.type = QMTYPE_SLIDER;	uiPlayerSetup.topColor.generic.flags = QMF_PULSEIFFOCUS|QMF_DROPSHADOW|addFlags;	uiPlayerSetup.topColor.generic.name = "Top color";	uiPlayerSetup.topColor.generic.x = 250;	uiPlayerSetup.topColor.generic.y = 550;	uiPlayerSetup.topColor.generic.width = 300;	uiPlayerSetup.topColor.generic.callback = UI_PlayerSetup_Callback;	uiPlayerSetup.topColor.generic.statusText = "Set a player model top color";	uiPlayerSetup.topColor.minValue = 0.0;	uiPlayerSetup.topColor.maxValue = 1.0;	uiPlayerSetup.topColor.range = 0.05f;	uiPlayerSetup.bottomColor.generic.id = ID_BOTTOMCOLOR;	uiPlayerSetup.bottomColor.generic.type = QMTYPE_SLIDER;	uiPlayerSetup.bottomColor.generic.flags = QMF_PULSEIFFOCUS|QMF_DROPSHADOW|addFlags;	uiPlayerSetup.bottomColor.generic.name = "Bottom color";	uiPlayerSetup.bottomColor.generic.x = 250;	uiPlayerSetup.bottomColor.generic.y = 620;	uiPlayerSetup.bottomColor.generic.width = 300;	uiPlayerSetup.bottomColor.generic.callback = UI_PlayerSetup_Callback;	uiPlayerSetup.bottomColor.generic.statusText = "Set a player model bottom color";	uiPlayerSetup.bottomColor.minValue = 0.0;	uiPlayerSetup.bottomColor.maxValue = 1.0;	uiPlayerSetup.bottomColor.range = 0.05f;	uiPlayerSetup.showModels.generic.id = ID_SHOWMODELS;	uiPlayerSetup.showModels.generic.type = QMTYPE_CHECKBOX;	uiPlayerSetup.showModels.generic.flags = QMF_HIGHLIGHTIFFOCUS|QMF_ACT_ONRELEASE|QMF_MOUSEONLY|QMF_DROPSHADOW|addFlags;	uiPlayerSetup.showModels.generic.name = "Show 3D Preview";	uiPlayerSetup.showModels.generic.x = 72;	uiPlayerSetup.showModels.generic.y = 380;	uiPlayerSetup.showModels.generic.callback = UI_PlayerSetup_Callback;	uiPlayerSetup.showModels.generic.statusText = "show 3D player models instead of preview thumbnails";	uiPlayerSetup.hiModels.generic.id = ID_HIMODELS;	uiPlayerSetup.hiModels.generic.type = QMTYPE_CHECKBOX;	uiPlayerSetup.hiModels.generic.flags = QMF_HIGHLIGHTIFFOCUS|QMF_ACT_ONRELEASE|QMF_MOUSEONLY|QMF_DROPSHADOW|addFlags;	uiPlayerSetup.hiModels.generic.name = "High quality models";	uiPlayerSetup.hiModels.generic.x = 72;	uiPlayerSetup.hiModels.generic.y = 430;	uiPlayerSetup.hiModels.generic.callback = UI_PlayerSetup_Callback;	uiPlayerSetup.hiModels.generic.statusText = "show hi-res models in multiplayer";	UI_PlayerSetup_GetConfig();	UI_AddItem( &uiPlayerSetup.menu, (void *)&uiPlayerSetup.background );	UI_AddItem( &uiPlayerSetup.menu, (void *)&uiPlayerSetup.banner );	UI_AddItem( &uiPlayerSetup.menu, (void *)&uiPlayerSetup.done );	UI_AddItem( &uiPlayerSetup.menu, (void *)&uiPlayerSetup.AdvOptions );	// disable playermodel preview for HLRally to prevent crash	if( game_hlRally == FALSE )		UI_AddItem( &uiPlayerSetup.menu, (void *)&uiPlayerSetup.view );	UI_AddItem( &uiPlayerSetup.menu, (void *)&uiPlayerSetup.name );	//if( !gMenu.m_gameinfo.flags & GFL_NOMODELS )	//{		UI_AddItem( &uiPlayerSetup.menu, (void *)&uiPlayerSetup.model );		UI_AddItem( &uiPlayerSetup.menu, (void *)&uiPlayerSetup.topColor );		UI_AddItem( &uiPlayerSetup.menu, (void *)&uiPlayerSetup.bottomColor );		UI_AddItem( &uiPlayerSetup.menu, (void *)&uiPlayerSetup.showModels );		UI_AddItem( &uiPlayerSetup.menu, (void *)&uiPlayerSetup.hiModels );	//}	// setup render and actor	uiPlayerSetup.refdef.fov_x = 40;	// NOTE: must be called after UI_AddItem whan we sure what UI_ScaleCoords is done	uiPlayerSetup.refdef.viewport[0] = uiPlayerSetup.view.generic.x;	uiPlayerSetup.refdef.viewport[1] = uiPlayerSetup.view.generic.y;	uiPlayerSetup.refdef.viewport[2] = uiPlayerSetup.view.generic.width;	uiPlayerSetup.refdef.viewport[3] = uiPlayerSetup.view.generic.height;	UI_PlayerSetup_CalcFov( &uiPlayerSetup.refdef );	uiPlayerSetup.ent = GET_MENU_EDICT ();	if( !uiPlayerSetup.ent )		return;	// adjust entity params	uiPlayerSetup.ent->curstate.number = 1;	// IMPORTANT: always set playerindex to 1	uiPlayerSetup.ent->curstate.animtime = gpGlobals->time;	// start animation	uiPlayerSetup.ent->curstate.sequence = 1;	uiPlayerSetup.ent->curstate.scale = 1.0f;	uiPlayerSetup.ent->curstate.frame = 0.0f;	uiPlayerSetup.ent->curstate.framerate = 1.0f;	uiPlayerSetup.ent->curstate.effects |= EF_LIGHT;	uiPlayerSetup.ent->curstate.controller[0] = 127;	uiPlayerSetup.ent->curstate.controller[1] = 127;	uiPlayerSetup.ent->curstate.controller[2] = 127;	uiPlayerSetup.ent->curstate.controller[3] = 127;	uiPlayerSetup.ent->latched.prevcontroller[0] = 127;	uiPlayerSetup.ent->latched.prevcontroller[1] = 127;	uiPlayerSetup.ent->latched.prevcontroller[2] = 127;	uiPlayerSetup.ent->latched.prevcontroller[3] = 127;	uiPlayerSetup.ent->origin[0] = uiPlayerSetup.ent->curstate.origin[0] = 45.0f / tan( DEG2RAD( uiPlayerSetup.refdef.fov_y / 2.0f ));	uiPlayerSetup.ent->origin[2] = uiPlayerSetup.ent->curstate.origin[2] = 2.0f;	uiPlayerSetup.ent->angles[1] = uiPlayerSetup.ent->curstate.angles[1] = 180.0f;	uiPlayerSetup.ent->player = true; // yes, draw me as playermodel}

bool TScriptInternalFunctions::runFunction(QString function, QStringList param, QString &result){    QString fn = function.toUpper();    if (fn == "ACOS") {        if (param.length() < 1)            return false;        bool ok = false;        double v = param[0].toDouble(&ok);        if (!ok)            v = 0;        result = QString::number(acos(v));        return true;    }    if (fn == "ASIN") {        if (param.length() < 1)            return false;        bool ok = false;        double v = param[0].toDouble(&ok);        if (!ok)            v = 0;        result = QString::number(asin(v));        return true;    }    if (fn == "ATAN") {        if (param.length() < 1)            return false;        bool ok = false;        double v = param[0].toDouble(&ok);        if (!ok)            v = 0;        result = QString::number(atan(v));        return true;    }    if (fn == "BUILDTYPE") {#ifdef STANDALONE        result = "STANDALONE";#endif#ifdef PACKAGED        result = "PACKAGED";#endif        return true;    }    if (fn == "CALC") {        // Calculate an expression (e.g. 5+5)        if (param.length() == 0)            return false;        QString expr = param[0];        result = calc(expr);        return true;    }    if (fn == "COS") {        if (param.length() < 1)            return false;        bool ok = false;        double v = param[0].toDouble(&ok);        if (!ok)            v = 0;        result = QString::number(cos(v));        return true;    }    if (fn == "COLORAT") { // $ColorAt(@window, layer, x, y)        if (param.length() < 3)            return false;        QString layer = "main";        if (param.length() > 3) {            layer = param[1];            param.removeAt(1);        }        subwindow_t sw = getCustomWindow(param[0]);        if (sw.type == WT_NOTHING)            return false;        int x = floor( param[1].toFloat() );        int y = floor( param[2].toFloat() );        result = sw.widget->picwinPtr()->colorAt(layer, x, y);        return true;    }    if (fn == "CURWINTYPE") {        // Returns the current target type (msg or channel)        subwindow_t sw = winList->value(*activeWid);//.........这里部分代码省略.........

double dTan(double x){  double v1 = tan(x);  return 1+v1*v1;}

/*===============UI_DrawPlayer===============*/void UI_DrawPlayer( float x, float y, float w, float h, uiPlayerInfo_t *pi, int time ) {	refdef_t		refdef;	refEntity_t		legs = {0};	refEntity_t		torso = {0};	refEntity_t		head = {0};	refEntity_t		gun = {0};	refEntity_t		barrel = {0};	refEntity_t		flash = {0};	vec3_t			origin;	int				renderfx;	vec3_t			mins = {-16, -16, -24};	vec3_t			maxs = {16, 16, 32};	float			len;	float			xx;	float			xscale;	float			yscale;	if ( !pi->legsModel || !pi->torsoModel || !pi->headModel || !pi->animations[0].numFrames ) {		return;	}	dp_realtime = time;	if ( pi->pendingWeapon != WP_NUM_WEAPONS && dp_realtime > pi->weaponTimer ) {		pi->weapon = pi->pendingWeapon;		pi->lastWeapon = pi->pendingWeapon;		pi->pendingWeapon = WP_NUM_WEAPONS;		pi->weaponTimer = 0;		if( pi->currentWeapon != pi->weapon ) {			trap_S_StartLocalSound( weaponChangeSound, CHAN_LOCAL );		}	}	CG_AdjustFrom640( &x, &y, &w, &h );	y -= jumpHeight;	memset( &refdef, 0, sizeof( refdef ) );	memset( &legs, 0, sizeof(legs) );	memset( &torso, 0, sizeof(torso) );	memset( &head, 0, sizeof(head) );	refdef.rdflags = RDF_NOWORLDMODEL;	AxisClear( refdef.viewaxis );	refdef.x = x;	refdef.y = y;	refdef.width = w;	refdef.height = h;	if ( ui_stretch.integer ) {		xscale = cgs.screenXScaleStretch;		yscale = cgs.screenYScaleStretch;	} else {		xscale = cgs.screenXScale;		yscale = cgs.screenYScale;	}	refdef.fov_x = (int)((float)refdef.width / xscale / 640.0f * 90.0f);	xx = refdef.width / xscale / tan( refdef.fov_x / 360 * M_PI );	refdef.fov_y = atan2( refdef.height / yscale, xx );	refdef.fov_y *= ( 360 / M_PI );	// calculate distance so the player nearly fills the box	len = 0.7 * ( maxs[2] - mins[2] );			origin[0] = len / tan( DEG2RAD(refdef.fov_x) * 0.5 );	origin[1] = 0.5 * ( mins[1] + maxs[1] );	origin[2] = -0.5 * ( mins[2] + maxs[2] );	refdef.time = dp_realtime;	trap_R_ClearScene();	// get the rotation information	UI_PlayerAngles( pi, legs.axis, torso.axis, head.axis );		// get the animation state (after rotation, to allow feet shuffle)	UI_PlayerAnimation( pi, &legs.oldframe, &legs.frame, &legs.backlerp,		 &torso.oldframe, &torso.frame, &torso.backlerp );	renderfx = RF_LIGHTING_ORIGIN | RF_NOSHADOW;	//	// add the legs	//	legs.hModel = pi->legsModel;	legs.customSkin = CG_AddSkinToFrame( &pi->modelSkin );	VectorCopy( origin, legs.origin );	VectorCopy( origin, legs.lightingOrigin );	legs.renderfx = renderfx;	VectorCopy (legs.origin, legs.oldorigin);//.........这里部分代码省略.........

double ddTan(double x){  double v1 = tan(x);  return 2*v1*(1+v1*v1);}

CAAPhysicalJupiterDetails CAAPhysicalJupiter::Calculate(double JD){  //What will be the return value  CAAPhysicalJupiterDetails details;  //Step 1  double d = JD - 2433282.5;  double T1 = d/36525;  double alpha0 = 268.00 + 0.1061*T1;  double alpha0rad = CAACoordinateTransformation::DegreesToRadians(alpha0);  double delta0 = 64.50 - 0.0164*T1;  double delta0rad = CAACoordinateTransformation::DegreesToRadians(delta0);  //Step 2  double W1 = CAACoordinateTransformation::MapTo0To360Range(17.710 + 877.90003539*d);  double W2 = CAACoordinateTransformation::MapTo0To360Range(16.838 + 870.27003539*d);    //Step 3  double l0 = CAAEarth::EclipticLongitude(JD);  double l0rad = CAACoordinateTransformation::DegreesToRadians(l0);  double b0 = CAAEarth::EclipticLatitude(JD);  double b0rad = CAACoordinateTransformation::DegreesToRadians(b0);  double R = CAAEarth::RadiusVector(JD);  //Step 4  double l = CAAJupiter::EclipticLongitude(JD);  double lrad = CAACoordinateTransformation::DegreesToRadians(l);  double b = CAAJupiter::EclipticLatitude(JD);  double brad = CAACoordinateTransformation::DegreesToRadians(b);  double r = CAAJupiter::RadiusVector(JD);  //Step 5  double x = r*cos(brad)*cos(lrad) - R*cos(l0rad);  double y = r*cos(brad)*sin(lrad) - R*sin(l0rad);  double z = r*sin(brad) - R*sin(b0rad);  double DELTA = sqrt(x*x + y*y + z*z);  //Step 6  l -= 0.012990*DELTA/(r*r);  lrad = CAACoordinateTransformation::DegreesToRadians(l);  //Step 7  x = r*cos(brad)*cos(lrad) - R*cos(l0rad);  y = r*cos(brad)*sin(lrad) - R*sin(l0rad);  z = r*sin(brad) - R*sin(b0rad);  DELTA = sqrt(x*x + y*y + z*z);  //Step 8  double e0 = CAANutation::MeanObliquityOfEcliptic(JD);  double e0rad = CAACoordinateTransformation::DegreesToRadians(e0);  //Step 9  double alphas = atan2(cos(e0rad)*sin(lrad) - sin(e0rad)*tan(brad), cos(lrad));  double deltas = asin(cos(e0rad)*sin(brad) + sin(e0rad)*cos(brad)*sin(lrad));  //Step 10  details.DS = CAACoordinateTransformation::RadiansToDegrees(asin(-sin(delta0rad)*sin(deltas) - cos(delta0rad)*cos(deltas)*cos(alpha0rad - alphas)));  //Step 11  double u = y*cos(e0rad) - z*sin(e0rad);  double v = y*sin(e0rad) + z*cos(e0rad);  double alpharad = atan2(u, x);  double alpha = CAACoordinateTransformation::RadiansToDegrees(alpharad);  double deltarad = atan2(v, sqrt(x*x + u*u));  double delta = CAACoordinateTransformation::RadiansToDegrees(deltarad);  double xi = atan2(sin(delta0rad)*cos(deltarad)*cos(alpha0rad - alpharad) - sin(deltarad)*cos(delta0rad), cos(deltarad)*sin(alpha0rad - alpharad));  //Step 12  details.DE = CAACoordinateTransformation::RadiansToDegrees(asin(-sin(delta0rad)*sin(deltarad) - cos(delta0rad)*cos(deltarad)*cos(alpha0rad - alpharad)));  //Step 13  details.Geometricw1 = CAACoordinateTransformation::MapTo0To360Range(W1 - CAACoordinateTransformation::RadiansToDegrees(xi) - 5.07033*DELTA);  details.Geometricw2 = CAACoordinateTransformation::MapTo0To360Range(W2 - CAACoordinateTransformation::RadiansToDegrees(xi) - 5.02626*DELTA);  //Step 14  double C = 57.2958 * (2*r*DELTA + R*R - r*r - DELTA*DELTA)/(4*r*DELTA);  if (sin(lrad - l0rad) > 0)  {    details.Apparentw1 = CAACoordinateTransformation::MapTo0To360Range(details.Geometricw1 + C);    details.Apparentw2 = CAACoordinateTransformation::MapTo0To360Range(details.Geometricw2 + C);  }  else  {    details.Apparentw1 = CAACoordinateTransformation::MapTo0To360Range(details.Geometricw1 - C);    details.Apparentw2 = CAACoordinateTransformation::MapTo0To360Range(details.Geometricw2 - C);  }  //Step 15  double NutationInLongitude = CAANutation::NutationInLongitude(JD);  double NutationInObliquity = CAANutation::NutationInObliquity(JD);  e0 += NutationInObliquity/3600;  e0rad = CAACoordinateTransformation::DegreesToRadians(e0);  //Step 16  alpha += 0.005693*(cos(alpharad)*cos(l0rad)*cos(e0rad) + sin(alpharad)*sin(l0rad))/cos(deltarad);  alpha = CAACoordinateTransformation::MapTo0To360Range(alpha);  alpharad = CAACoordinateTransformation::DegreesToRadians(alpha);  delta += 0.005693*(cos(l0rad)*cos(e0rad)*(tan(e0rad)*cos(deltarad) - sin(alpharad)*sin(deltarad)) + cos(alpharad)*sin(deltarad)*sin(l0rad));  deltarad = CAACoordinateTransformation::DegreesToRadians(delta);//.........这里部分代码省略.........

void FollowCamera::update(Step * _step){		lastOrientation = childTransform->getOrientationQuat();	glm::quat newOrientation = glm::quat(1.f, 0.f, 0.f, 0.f);	newOrientation = glm::rotate(newOrientation, yaw, upVectorLocal);	newOrientation = glm::rotate(newOrientation, pitch, rightVectorLocal);	newOrientation = glm::slerp(lastOrientation, newOrientation, 0.15f * static_cast<float>(sweet::deltaTimeCorrection));	childTransform->setOrientation(newOrientation);	forwardVectorRotated   = newOrientation * forwardVectorLocal;	rightVectorRotated	   = newOrientation * rightVectorLocal;	upVectorRotated		   = newOrientation * upVectorLocal;		lookAtSpot = glm::vec3(0.f, 0.f, 0.f);	float targetMinX = 9999999999.f;	float targetMinY = 9999999999.f;	float targetMaxX = -9999999999.f;	float targetMaxY = -9999999999.f;		for(signed long int i = targets.size()-1; i >= 0; --i){		if(!targets.at(i).active){			if(targets.at(i).weight <= 0.001f){				targets.erase(targets.begin() + i);			}		}else{			targets.at(i).pos = targets.at(i).target->getWorldPos();		}	}	for(Target & t : targets){		targetMinX = std::min((t.pos.x-buffer)*t.weight, targetMinX);		targetMaxX = std::max((t.pos.x+buffer)*t.weight, targetMaxX);		targetMinY = std::min((t.pos.y-buffer)*t.weight, targetMinY);		targetMaxY = std::max((t.pos.y+buffer)*t.weight, targetMaxY);		if(t.active){			t.weight = std::min(1.f, t.weight + 0.05f);		}else{			t.weight = std::max(0.f, t.weight - 0.01f);		}	}	float screenWidth = targetMaxX - targetMinX;	float screenHeight = targetMaxY - targetMinY;		// move camera	lookAtSpot.x = targetMinX;	lookAtSpot.y = targetMinY;		lookAtSpot += offset;		if(useBounds){		if(minBounds.height != 0){			if(lookAtSpot.y < minBounds.y){				lookAtSpot.y = minBounds.y;			}			if(lookAtSpot.y + screenHeight > minBounds.x + minBounds.height){				lookAtSpot.y -= (lookAtSpot.y + screenHeight - (minBounds.y + minBounds.height));			}			if(lookAtSpot.y < minBounds.y){				screenHeight -= minBounds.y - lookAtSpot.y;				lookAtSpot.y = minBounds.y;			}		}		if(minBounds.width != 0){			if(lookAtSpot.x < minBounds.x){				lookAtSpot.x = minBounds.x;			}			if(lookAtSpot.x + screenWidth > minBounds.x + minBounds.width){				lookAtSpot.x -= (lookAtSpot.x + screenWidth - (minBounds.x + minBounds.width));			}			if(lookAtSpot.x < minBounds.x){				screenWidth -= minBounds.x - lookAtSpot.x;				lookAtSpot.x = minBounds.x;			}		}	}	// calculate zoom and account for FoV (the camera FoV seems to be vertical, so if the screen w > screen h, we need to take the h / the intended aspect ratio)	float ar1 = screenWidth/screenHeight;	glm::vec2 screenDimensions = sweet::getWindowDimensions();	float ar2 = static_cast<float>(screenDimensions.x)/static_cast<float>(screenDimensions.y);	float zoom;	if(ar1 > ar2){		zoom = std::max(minimumZoom, screenWidth / ar2);	}else if(ar1 < ar2){		zoom = std::max(minimumZoom, screenHeight);	}else{		zoom = std::max(minimumZoom, screenHeight);	}	lookAtSpot.x += screenWidth * 0.5f;	lookAtSpot.y += screenHeight * 0.5f;	lookAtSpot += offset;	float dist = zoom / (tan(glm::radians(fieldOfView) * 0.5f) * 2.f);//.........这里部分代码省略.........

float Light::MicroFacet(Vector l, Vector v, Vector n, float m) {    Vector h = (l + v) / Norm(l + v);    float J = acos(DotProduct(h, n));    auto Beckmann = static_cast<float>((1 / (4.0f * (m * m) * pow((cos(J)), 4))) * exp((-1.0f * (tan(J) * tan(J))) / (m * m)));    return Beckmann;}

void latlon2_(double *alat1, double *alon1, double *delta, double *azi, double *alat2, double *alon2) {	double alat, alatr, alon, b, c, coslat, dlon;	double r13, sinlat, x1, x2, x3;	/* changed for ellipticity of earth	 * changed use of *alat1 and *alat2	 */	double esq, alat3;	esq=(1.0-1.0/298.25)*(1.0-1.0/298.25);	alat3=atan(tan(*alat1*DEG_TO_RAD)*esq)*RAD_TO_DEG;	/* Convert a geographical location to geocentric cartesian 	 * coordinates, assuming a spherical earth	 */	alat = 90.0 - *delta;	alon = 180.0 - *azi;	r13  = cos(DEG_TO_RAD*alat);	/* x1:	Axis 1 intersects equator at  0 deg longitude  	 * x2:	Axis 2 intersects equator at 90 deg longitude  	 * x3:	Axis 3 intersects north pole	 */	x1 = r13*sin(DEG_TO_RAD*alon);	x2 = sin(DEG_TO_RAD*alat);	x3 = r13*cos(DEG_TO_RAD*alon);	/* Rotate in cartesian coordinates.  The cartesian coordinate system 	 * is most easily described in geographic terms.  The origin is at 	 * the Earth's center.  Rotation by alat1 degrees southward, about 	 * the 1-axis.	 */	alatr  = (90.0-alat3)/RAD_TO_DEG;	sinlat = sin(alatr);	coslat = cos(alatr);	b      = x2;	c      = x3;	x2     = b*coslat - c*sinlat;	x3     = b*sinlat + c*coslat;	/* Convert geocentric cartesian coordinates to a geographical 	 * location, assuming a spherical earth	 */	r13    = sqrt(x3*x3 + x1*x1);	dlon   = RAD_TO_DEG*atan2(x1, x3);	/*  changed for ellipticity of earth	 *   *alat2 = RAD_TO_DEG*atan2(x2, r13);	 */	alat3= atan2(x2, r13);	*alat2=RAD_TO_DEG * atan(tan(alat3)/esq);	*alon2 = *alon1 + dlon;	if (fabs(*alon2) > 180.0)		*alon2 = SIGN((360.0-fabs(*alon2)), *alon2);}

Catoms2DSimulator::Catoms2DSimulator(int argc, char *argv[], Catoms2DBlockCode *(*catoms2DBlockCodeBuildingFunction)(Catoms2DBlock*)) : BaseSimulator::Simulator(argc, argv) {	OUTPUT << "/033[1;34m" << "Catoms2DSimulator constructor" << "/033[0m" << endl;	int currentID = 1;	Catoms2DWorld *world = NULL;	buildNewBlockCode = catoms2DBlockCodeBuildingFunction;	float blockSize[3];	testMode = false;	/* reading the xml file */	TiXmlNode *node = xmlDoc->FirstChild("world");	if (node) {		TiXmlElement* worldElement = node->ToElement();		const char *attr= worldElement->Attribute("gridSize");		int lx,ly,lz;		if (attr) {			string str=attr;			int pos = str.find_first_of(',');			lx = atoi(str.substr(0,pos).c_str());			ly = 1;			lz = atoi(str.substr(pos+1,str.length()-pos-1).c_str());			OUTPUT << "grid size : " << lx << " x " << ly << " x " << lz << endl;		} else {			OUTPUT << "WARNING No grid size in XML file" << endl;		}		attr=worldElement->Attribute("windowSize");		if (attr) {			string str=attr;	 		int pos = str.find_first_of(',');			GlutContext::initialScreenWidth = atoi(str.substr(0,pos).c_str());			GlutContext::initialScreenHeight = atoi(str.substr(pos+1,str.length()-pos-1).c_str());			GlutContext::screenWidth = GlutContext::initialScreenWidth;			GlutContext::screenHeight = GlutContext::initialScreenHeight;		}		createWorld(lx, ly, lz, argc, argv);		world = getWorld();		world->loadTextures("../../simulatorCore/catoms2DTextures");	} else {		ERRPUT << "ERROR : NO world in XML file" << endl;		exit(1);	}	createScheduler();	// loading the camera parameters	TiXmlNode *nodeConfig = node->FirstChild("camera");	if (nodeConfig) {		TiXmlElement* cameraElement = nodeConfig->ToElement();		const char *attr=cameraElement->Attribute("target");		double def_near=1,def_far=1500;		float angle=45.0;		if (attr) {			string str(attr);			int pos = str.find_first_of(',');			Vecteur target;			target.pt[0] = atof(str.substr(0,pos).c_str());			target.pt[1] = 1;			target.pt[2] = atoi(str.substr(pos+1,str.length()-pos-1).c_str());			world->getCamera()->setTarget(target);		}		attr=cameraElement->Attribute("angle");		if (attr) {			angle = atof(attr);			world->getCamera()->setAngle(angle);		}		attr=cameraElement->Attribute("directionSpherical");		if (attr) {			string str(attr);			int pos1 = str.find_first_of(','),			pos2 = str.find_last_of(',');			float az,ele,dist;			az = -90.0+atof(str.substr(0,pos1).c_str());			ele = atof(str.substr(pos1+1,pos2-pos1-1).c_str());			dist = atof(str.substr(pos2+1,str.length()-pos1-1).c_str());			world->getCamera()->setDirection(az,ele);			world->getCamera()->setDistance(dist);			az = dist*sin(angle*M_PI/180.0);			def_near = dist-az;			def_far = dist+az;		}		attr=cameraElement->Attribute("near");		if (attr) {			def_near = atof(attr);		}		attr=cameraElement->Attribute("far");		if (attr) {			def_far = atof(attr);		}		world->getCamera()->setNearFar(def_near,def_far);	}	// loading the spotlight parameters	nodeConfig = node->FirstChild("spotlight");	if (nodeConfig) {		Vecteur target;		float az=0,ele=60,dist=1000,angle=50;		TiXmlElement* lightElement = nodeConfig->ToElement();//.........这里部分代码省略.........

double f( double x ){      return ( p * exp( -x ) + q*sin( x ) + r*cos( x ) + s*tan( x ) + t*x*x + u );}

bool GeoAlgorithms::initVincenty(double aLat1, double aLon1, double aLat2, double aLon2){   // Verify that input latitudes are between -90 and 90 and longitudes are   // between -180 and 180   if ((abs(aLat1) > 90) || (abs(aLat2) > 90) || (abs(aLon1) > 180)        || (abs(aLon2) > 180))   {       return false;   }   // convert inputs in degrees to radians:   aLat1 = aLat1 * 0.0174532925199433;   aLon1 = aLon1 * 0.0174532925199433;   aLat2 = aLat2 * 0.0174532925199433;   aLon2 = aLon2 * 0.0174532925199433;   // correct for errors at exact poles by adjusting 0.6 millimeters:   if (abs(GeoConversions::PI_OVER_2-abs(aLat1)) < (1e-10))   {       aLat1 = getSign(aLat1) * (GeoConversions::PI_OVER_2 - (1e-10));   }   if (abs(GeoConversions::PI_OVER_2-abs(aLat2)) < (1e-10))   {       aLat2 = getSign(aLat2) * (GeoConversions::PI_OVER_2 - (1e-10));   }   // Ellipse CalcuaAltitudeions?   mVincentyU1 = atan(m1MinF*tan(aLat1));   mVincentyU2 = atan(m1MinF*tan(aLat2));   aLon1 = getMod(aLon1, (GeoConversions::TWO_PI));   aLon2 = getMod(aLon2, (GeoConversions::TWO_PI));   mVincentyL = aLon2-aLon1;   if (abs(mVincentyL) > PI)   {     mVincentyL = getSign(mVincentyL) * (GeoConversions::TWO_PI - abs(mVincentyL));   }   // Initialize Variables for Loop   double sin_mVincentyU1 = sin(mVincentyU1);   double cos_mVincentyU1 = cos(mVincentyU1);   double sin_mVincentyU2 = sin(mVincentyU2);   double cos_mVincentyU2 = cos(mVincentyU2);   double sinU1_sinU2 = sin_mVincentyU1 * sin_mVincentyU2;   double cosU1_sinU2 = cos_mVincentyU1 * sin_mVincentyU2;   double sinU1_cosU2 = sin_mVincentyU1 * cos_mVincentyU2;   double cosU1_cosU2 = cos_mVincentyU1 * cos_mVincentyU2;   double sin_mVincentyLambda = 0;   double cos_mVincentyLambda = 0;   double cos_mVincentyAlpha = 0;   double sin_mVincentySigma = 0;   double cos_mVincentySigma = 0;   double lLambdaOld = 0;   long   lIterCount = 0;   double lSinSigma = 0;   double lCosSigma = 0;   double c = 0;   mVincentySigma = 0;   mVincentyAlpha = 0;   mVincentyCosToSigmaM = 0;   mVincentyLambda = mVincentyL;   // ?   while ((!lIterCount) || abs((mVincentyLambda-lLambdaOld) > (1e-12)))   {      lIterCount += 1;      if (lIterCount > 50)      {          mVincentyLambda = PI;          break;      }      sin_mVincentyLambda = sin(mVincentyLambda);      cos_mVincentyLambda = cos(mVincentyLambda);      lLambdaOld = mVincentyLambda;      lSinSigma = sqrt(pow(cos_mVincentyU2 * sin_mVincentyLambda, 2) +         pow(cosU1_sinU2 - sinU1_cosU2 * cos_mVincentyLambda, 2));      lCosSigma = sinU1_sinU2 + cosU1_cosU2 * cos_mVincentyLambda;      mVincentySigma = atan2(lSinSigma, lCosSigma);      sin_mVincentySigma = sin(mVincentySigma);      cos_mVincentySigma = cos(mVincentySigma);      mVincentyAlpha = asin(cosU1_cosU2 * sin_mVincentyLambda / sin_mVincentySigma);      cos_mVincentyAlpha = cos(mVincentyAlpha);      mVincentyCosToSigmaM = cos_mVincentySigma -  2 * sinU1_sinU2 / pow(cos_mVincentyAlpha, 2);      c = mF/ 16 * pow(cos_mVincentyAlpha, 2) * (4 + mF * (4 - 3 * pow(cos_mVincentyAlpha, 2)));      mVincentyLambda = mVincentyL + (1 - c) * mF * sin(mVincentyAlpha) * (mVincentySigma + c * sin_mVincentySigma *         (mVincentyCosToSigmaM + c * cos_mVincentySigma * (-1 + 2 * pow(mVincentyCosToSigmaM, 2))));      // Correct for convergence failure in the case of essentially      // antipodal points      if (mVincentyLambda > PI)      {          mVincentyLambda = PI;          break;      }   }   return true;}

//---------------------int fdct_wrapping_invsepangle(double XL1, double XL2, int nbangle, vector<CpxNumMat>& csc, CpxOffMat& Xhgh){  typedef pair<int,int> intpair;  map<intpair, fftwnd_plan> planmap;    int XS1, XS2;  int XF1, XF2;  double XR1, XR2;	 fdct_wrapping_rangecompute(XL1, XL2, XS1, XS2, XF1, XF2, XR1, XR2);  Xhgh.resize(XS1, XS2);    int nbquadrants = 4;  int nd = nbangle / 4;  int wcnt = 0;    //backup  CpxOffMat Xhghb(Xhgh);  double XL1b = XL1;  double XL2b = XL2;  int qvec[] = {2,1,0,3};  for(int qi=0; qi<nbquadrants; qi++) {	 int q = qvec[qi];	 //ROTATE data to its right position	 fdct_wrapping_rotate_forward(q, XL1b, XL2b, XL1, XL2);	 XL1 = abs(XL1);	 XL2 = abs(XL2);	 fdct_wrapping_rotate_forward(q, Xhghb, Xhgh);	 //figure out XS, XF, XR	 double XW1 = XL1/nd;	 double XW2 = XL2/nd;	 int XS1, XS2;  int XF1, XF2;  double XR1, XR2;  fdct_wrapping_rangecompute(XL1, XL2, XS1, XS2, XF1, XF2, XR1, XR2);	 for(int w=nd-1; w>=0; w--) {		double xs = XR1/4 - (XW1/2)/4;		double xe = XR1;		double ys = -XR2 + (w-0.5)*XW2;		double ye = -XR2 + (w+1.5)*XW2; //x range		int xn = int(ceil(xe-xs));			 int yn = int(ceil(ye-ys));		//MAKE THEM ODD		if(xn%2==0) xn++;		if(yn%2==0) yn++;		int xf = int(ceil(xs));			 //int yf = int(ceil(ys));		//theta		double thts, thtm, thte; //y direction		if(w==0) {		  thts = atan2(-1.0, 1.0-1.0/nd);		  thtm = atan2(-1.0+1.0/nd, 1.0);		  thte = atan2(-1.0+3.0/nd, 1.0);		} else if(w==nd-1) {		  thts = atan2(-1.0+(2.0*w-1.0)/nd, 1.0);		  thtm = atan2(-1.0+(2.0*w+1.0)/nd, 1.0);		  thte = atan2(1.0, 1.0-1.0/nd);		} else {		  thts = atan2(-1.0+(2.0*w-1.0)/nd, 1.0);		  thtm = atan2(-1.0+(2.0*w+1.0)/nd, 1.0);		  thte = atan2(-1.0+(2.0*w+3.0)/nd, 1.0);		}		int xh = xn/2;		int yh = yn/2; //half length		CpxOffMat wpdata(xn,yn);		{		  //load		  int xn = csc[wcnt].m();		  int yn = csc[wcnt].n();		  CpxNumMat tpdata(csc[wcnt]);		  //fft		  fftwnd_plan p = NULL;		  map<intpair, fftwnd_plan>::iterator mit=planmap.find( intpair(xn,yn) );		  if(mit!=planmap.end()) {			 p = (*mit).second;		  } else {			 p = fftw2d_create_plan(yn, xn, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_IN_PLACE);			 planmap[ intpair(xn, yn) ] = p;		  }		  fftwnd_one(p, (fftw_complex*)tpdata.data(), NULL);		  double sqrtprod = sqrt(double(xn*yn));		  for(int i=0; i<xn; i++)		  for(int j=0; j<yn; j++)			 tpdata(i,j) /= sqrtprod;		  //fftshift		  CpxOffMat rpdata;		  fdct_wrapping_fftshift(tpdata,rpdata);		  //rotate forward		  fdct_wrapping_rotate_forward(q, rpdata, wpdata);		}				double R21 = XR2/XR1; //ratio		for(int xcur=xf; xcur<xe; xcur++) { //for each layer		  int yfm = (int)ceil( max(-XR2, R21*xcur*tan(thts)) );		  int yto = (int)floor( min(XR2, R21*xcur*tan(thte)) );		  for(int ycur=yfm; ycur<=yto; ycur++) {			 int tmpx = xcur%xn;				  if(tmpx<-xh) tmpx+=xn;				  if(tmpx>=-xh+xn) tmpx-=xn;			 int tmpy = ycur%yn;				  if(tmpy<-yh) tmpy+=yn;				  if(tmpy>=-yh+yn) tmpy-=yn;			 //partition of unity			 double thtcur = atan2(ycur/XR2, xcur/XR1);			 double wtht;			 if(thtcur<thtm) {				double l,r; fdct_wrapping_window((thtcur-thts)/(thtm-thts), l, r);				wtht = l;			 } else {				double l,r; fdct_wrapping_window((thtcur-thtm)/(thte-thtm), l, r);				wtht = r;			 }			 double pou = wtht;			 wpdata(tmpx,tmpy) *= pou;			 Xhgh(xcur,ycur) += wpdata(tmpx,tmpy);		  }		}				wcnt++;//.........这里部分代码省略.........

//.........这里部分代码省略.........					trap_R_SetFog( FOG_PORTALVIEW, fogStart, fogEnd, fogColor[0], fogColor[1], fogColor[2], 1.1 );					foginited = qtrue;				}			} else {				if ( !foginited ) {					trap_R_SetFog( FOG_PORTALVIEW, 0,0,0,0,0,0 ); // init to null					foginited = qtrue;				}			}		}		//----(SA)	end		if ( cg.predictedPlayerState.pm_type == PM_INTERMISSION ) {			// if in intermission, use a fixed value			fov_x = 90;		} else {			// user selectable			if ( cgs.dmflags & DF_FIXED_FOV ) {				// dmflag to prevent wide fov for all clients				fov_x = 90;			} else {				fov_x = cg_fov.value;				if ( fov_x < 1 ) {					fov_x = 1;				} else if ( fov_x > 160 ) {					fov_x = 160;				}			}			// account for zooms			if ( cg.zoomval ) {				zoomFov = cg.zoomval;   // (SA) use user scrolled amount				if ( zoomFov < 1 ) {					zoomFov = 1;				} else if ( zoomFov > 160 ) {					zoomFov = 160;				}			} else {				zoomFov = lastfov;			}			// do smooth transitions for the binocs			if ( cg.zoomedBinoc ) {        // binoc zooming in				f = ( cg.time - cg.zoomTime ) / (float)ZOOM_TIME;				if ( f > 1.0 ) {					fov_x = zoomFov;				} else {					fov_x = fov_x + f * ( zoomFov - fov_x );				}				lastfov = fov_x;			} else if ( cg.zoomval ) {    // zoomed by sniper/snooper				fov_x = cg.zoomval;				lastfov = fov_x;			} else {                    // binoc zooming out				f = ( cg.time - cg.zoomTime ) / (float)ZOOM_TIME;				if ( f > 1.0 ) {					fov_x = fov_x;				} else {					fov_x = zoomFov + f * ( fov_x - zoomFov );				}			}		}		if ( cg.weaponSelect == WP_SNOOPERSCOPE ) {			cg.refdef.rdflags |= RDF_SNOOPERVIEW;		} else {			cg.refdef.rdflags &= ~RDF_SNOOPERVIEW;		}		if ( cg.snap->ps.persistant[PERS_HWEAPON_USE] ) {			fov_x = 55;		}		x = cg.refdef.width / tan( fov_x / 360 * M_PI );		fov_y = atan2( cg.refdef.height, x );		fov_y = fov_y * 360 / M_PI;		cg.refdef.fov_x = fov_x;		cg.refdef.fov_y = fov_y;		cg.refdef.rdflags |= RDF_SKYBOXPORTAL;		cg.refdef.rdflags |= RDF_DRAWSKYBOX;	} else {    // end if(cg_skybox.integer)		cg.refdef.rdflags |= RDF_SKYBOXPORTAL;		cg.refdef.rdflags &= ~RDF_DRAWSKYBOX;	}	cg.refdef.time = cg.time;	// draw the skybox	trap_R_RenderScene( &cg.refdef );	cg.refdef = backuprefdef;}

void Game::setupViewpoint(SIDE side){	//22.5 correspond à l'angle de vision du viewport divisé par 2, en degrés	setupViewpoint(		0.0,		0.0,		static_cast<float>(side * (world->getDepth() / 2 + Tools<float>::maximum(world->getHeight(), world->getWidth()) / (2 * tan(22.5 * Tools<int>::pi() / 180)))),		0.0,		0.0,		1.0f		);}

void turret_breach_think(edict_t * self){    edict_t *ent;    vec3_t current_angles;    vec3_t delta;    VectorCopy(self->s.angles, current_angles);    AnglesNormalize(current_angles);    AnglesNormalize(self->move_angles);    if (self->move_angles[PITCH] > 180)        self->move_angles[PITCH] -= 360;    // clamp angles to mins & maxs    if (self->move_angles[PITCH] > self->pos1[PITCH])        self->move_angles[PITCH] = self->pos1[PITCH];    else if (self->move_angles[PITCH] < self->pos2[PITCH])        self->move_angles[PITCH] = self->pos2[PITCH];    if ((self->move_angles[YAW] < self->pos1[YAW])            || (self->move_angles[YAW] > self->pos2[YAW])) {        float dmin, dmax;        dmin = fabs(self->pos1[YAW] - self->move_angles[YAW]);        if (dmin < -180)            dmin += 360;        else if (dmin > 180)            dmin -= 360;        dmax = fabs(self->pos2[YAW] - self->move_angles[YAW]);        if (dmax < -180)            dmax += 360;        else if (dmax > 180)            dmax -= 360;        if (fabs(dmin) < fabs(dmax))            self->move_angles[YAW] = self->pos1[YAW];        else            self->move_angles[YAW] = self->pos2[YAW];    }    VectorSubtract(self->move_angles, current_angles, delta);    if (delta[0] < -180)        delta[0] += 360;    else if (delta[0] > 180)        delta[0] -= 360;    if (delta[1] < -180)        delta[1] += 360;    else if (delta[1] > 180)        delta[1] -= 360;    delta[2] = 0;    if (delta[0] > self->speed * FRAMETIME)        delta[0] = self->speed * FRAMETIME;    if (delta[0] < -1 * self->speed * FRAMETIME)        delta[0] = -1 * self->speed * FRAMETIME;    if (delta[1] > self->speed * FRAMETIME)        delta[1] = self->speed * FRAMETIME;    if (delta[1] < -1 * self->speed * FRAMETIME)        delta[1] = -1 * self->speed * FRAMETIME;    VectorScale(delta, 1.0 / FRAMETIME, self->avelocity);    self->nextthink = level.time + FRAMETIME;    for (ent = self->teammaster; ent; ent = ent->teamchain)        ent->avelocity[1] = self->avelocity[1];    // if we have adriver, adjust his velocities    if (self->owner) {        float angle;        float target_z;        float diff;        vec3_t target;        vec3_t dir;        // angular is easy, just copy ours        self->owner->avelocity[0] = self->avelocity[0];        self->owner->avelocity[1] = self->avelocity[1];        // x & y        angle = self->s.angles[1] + self->owner->move_origin[1];        angle *= (M_PI * 2 / 360);        target[0] =            SnapToEights(self->s.origin[0] +                         cos(angle) * self->owner->move_origin[0]);        target[1] =            SnapToEights(self->s.origin[1] +                         sin(angle) * self->owner->move_origin[0]);        target[2] = self->owner->s.origin[2];        VectorSubtract(target, self->owner->s.origin, dir);        self->owner->velocity[0] = dir[0] * 1.0 / FRAMETIME;        self->owner->velocity[1] = dir[1] * 1.0 / FRAMETIME;        // z        angle = self->s.angles[PITCH] * (M_PI * 2 / 360);        target_z =            SnapToEights(self->s.origin[2] +                         self->owner->move_origin[0] * tan(angle) +                         self->owner->move_origin[2]);//.........这里部分代码省略.........

/** * /brief Calculate distance between two points * This function uses an algorithm for an oblate spheroid earth model. * The algorithm is described here:  * http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf * /return Distance in meters */double nmea_distance_ellipsoid(        const nmeaPOS *from_pos,    /**< From position in radians */        const nmeaPOS *to_pos,      /**< To position in radians */        double *from_azimuth,       /**< (O) azimuth at "from" position in radians */        double *to_azimuth          /**< (O) azimuth at "to" position in radians */        ){    /* All variables */    double f, a, b, sqr_a, sqr_b;    double L, phi1, phi2, U1, U2, sin_U1, sin_U2, cos_U1, cos_U2;    double sigma, sin_sigma, cos_sigma, cos_2_sigmam, sqr_cos_2_sigmam, sqr_cos_alpha, lambda, sin_lambda, cos_lambda, delta_lambda;    int remaining_steps;     double sqr_u, A, B, delta_sigma;    /* Check input */    NMEA_ASSERT(from_pos != 0);    NMEA_ASSERT(to_pos != 0);    if ((from_pos->lat == to_pos->lat) && (from_pos->lon == to_pos->lon))    { /* Identical points */        if ( from_azimuth != 0 )            *from_azimuth = 0;        if ( to_azimuth != 0 )            *to_azimuth = 0;        return 0;        } /* Identical points */    /* Earth geometry */    f = NMEA_EARTH_FLATTENING;    a = NMEA_EARTH_SEMIMAJORAXIS_M;    b = (1 - f) * a;    sqr_a = a * a;    sqr_b = b * b;    /* Calculation */    L = to_pos->lon - from_pos->lon;    phi1 = from_pos->lat;    phi2 = to_pos->lat;    U1 = atan((1 - f) * tan(phi1));    U2 = atan((1 - f) * tan(phi2));    sin_U1 = sin(U1);    sin_U2 = sin(U2);    cos_U1 = cos(U1);    cos_U2 = cos(U2);    /* Initialize iteration */    sigma = 0;    sin_sigma = sin(sigma);    cos_sigma = cos(sigma);    cos_2_sigmam = 0;    sqr_cos_2_sigmam = cos_2_sigmam * cos_2_sigmam;    sqr_cos_alpha = 0;    lambda = L;    sin_lambda = sin(lambda);                                cos_lambda = cos(lambda);                           delta_lambda = lambda;    remaining_steps = 20;     while ((delta_lambda > 1e-12) && (remaining_steps > 0))     { /* Iterate */        /* Variables */        double tmp1, tmp2, sin_alpha, cos_alpha, C, lambda_prev;        /* Calculation */        tmp1 = cos_U2 * sin_lambda;        tmp2 = cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda;          sin_sigma = sqrt(tmp1 * tmp1 + tmp2 * tmp2);                        cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda;           sin_alpha = cos_U1 * cos_U2 * sin_lambda / sin_sigma;          cos_alpha = cos(asin(sin_alpha));                         sqr_cos_alpha = cos_alpha * cos_alpha;                             cos_2_sigmam = cos_sigma - 2 * sin_U1 * sin_U2 / sqr_cos_alpha;        sqr_cos_2_sigmam = cos_2_sigmam * cos_2_sigmam;         C = f / 16 * sqr_cos_alpha * (4 + f * (4 - 3 * sqr_cos_alpha));        lambda_prev = lambda;         sigma = asin(sin_sigma);         lambda = L +             (1 - C) * f * sin_alpha            * (sigma + C * sin_sigma * (cos_2_sigmam + C * cos_sigma * (-1 + 2 * sqr_cos_2_sigmam)));                                                        delta_lambda = lambda_prev - lambda;         if ( delta_lambda < 0 ) delta_lambda = -delta_lambda;         sin_lambda = sin(lambda);        cos_lambda = cos(lambda);        remaining_steps--;     }  /* Iterate */    /* More calculation  */    sqr_u = sqr_cos_alpha * (sqr_a - sqr_b) / sqr_b;     A = 1 + sqr_u / 16384 * (4096 + sqr_u * (-768 + sqr_u * (320 - 175 * sqr_u)));    B = sqr_u / 1024 * (256 + sqr_u * (-128 + sqr_u * (74 - 47 * sqr_u)));    delta_sigma = B * sin_sigma * (         cos_2_sigmam + B / 4 * (         cos_sigma * (-1 + 2 * sqr_cos_2_sigmam) -//.........这里部分代码省略.........