C++ cosf函数代码示例

static inline void calc_coeffs4 (double *bcoeff, double *acoeff, double w0, double Dw, double gain)
{
  double G = exp10f(gain/20.f);
  // bandwidth gain is set to half gain
  double GB = exp10f(gain/40.f);
  double c0 = cosf(w0);
  double WB = tanf(Dw/2);
  double e = sqrtf((G*G-GB*GB)/(GB*GB-1.f));
  double g = powf(G, 0.25f);
  double a = powf(e, 0.25f);
#define  si1  0.382683432365089781779f;
#define  si2  0.923879532511286738483f;
  //Ba(1+i,:) = [g^2*WB^2*v, 2*g*b*si*WB, b^2*v]
  //Aa(1+i,:) = [WB^2*v, 2*a*si*WB, a^2*v]
  double Aa0 = WB*WB;
  double Aatmp = 2*a*WB;
  double Aa10 = Aatmp * si1;
  double Aa11 = Aatmp * si2;
  double as = a*a;
  double Ba0 = Aa0 * g*g;
  double Batmp = 2*g * a * WB;
  double Ba10 = Batmp * si1;
  double Ba11 = Batmp * si2;
  //  fprintf(stdout, "%f %f %f %f %f %f %f\n", Aa0, Aa10, Aa11, Ba0, Ba10, Ba11, as);
  //  double Bhat0 = Ba0/Aa0;
  //  D = A0(i)+A1(i)+A2(i)
  double D1 = Aa0+Aa11+as;
  double D0 = Aa0+Aa10+as;
  //Bhat(i,1) = (B0(i)+B1(i)+B2(i))./D
  //(B0(i)+B1(i)+B2(i))./D
  double tmp = Ba0+as;
  double x;
  double Bhat00 = (Ba10+tmp)/D0;
  double Bhat01 = (Ba11+tmp)/D1;
  x=2.f*(Ba0-as);
  double Bhat10 = x/D0;
  double Bhat11 = x/D1;
  double Bhat20 = (tmp-Ba10)/D0;
  double Bhat21 = (tmp-Ba11)/D1;
  //Ahat(i,1) = 1
  //Ahat(i,2) = 2*(A0(i)-A2(i))./D
  //Ahat(i,3) = (A0(i)-A1(i)+A2(i))./D
  x = 2.f*(Aa0-as);
  double Ahat10 = x/D0;
  double Ahat11 = x/D1;
  double Ahat20 = (Aa0-Aa10+as)/D0;
  double Ahat21 = (Aa0-Aa11+as)/D1;

  //  fprintf(stdout, "%f %f %f %f %f\n %f %f %f %f\n\n",Bhat00, Bhat01, Bhat10, Bhat11, Bhat21, Ahat11, Ahat10, Ahat20, Ahat21);


  // B(i,2) = c0*(Bhat(i,2)-2*Bhat(i,1))
  bcoeff[0] = Bhat00;
  bcoeff[5] = Bhat01;
  bcoeff[1] = c0*(Bhat10-2.f*Bhat00);
  bcoeff[6] = c0*(Bhat11-2.f*Bhat01);
  bcoeff[2] = (Bhat00-Bhat10+Bhat20)*c0*c0- Bhat10;
  bcoeff[7] = (Bhat01-Bhat11+Bhat21)*c0*c0 - Bhat11;
  bcoeff[3] = c0*(Bhat10-2.f*Bhat20);
  bcoeff[8] = c0*(Bhat11-2.f*Bhat21);
  bcoeff[4] = Bhat20;
  bcoeff[9] = Bhat21;

  //  fprintf(stdout, "%f %f %f %f %f\n %f %f %f %f %f\n", Bres00, Bres10, Bres20, Bres30, Bres40, Bres01, Bres11, Bres21, Bres31, Bres41);

  acoeff[0]= c0*(Ahat10-2.f);
  acoeff[4]= c0*(Ahat11-2.f);
  acoeff[1]= (1.f-Ahat10+Ahat20)*c0*c0 - Ahat10;
  acoeff[5]= (1.f-Ahat11+Ahat21)*c0*c0- Ahat11;
  acoeff[2]= c0*(Ahat10-2.f*Ahat20);
  acoeff[6]= c0*(Ahat11-2.f*Ahat21);
  acoeff[3]= Ahat20;
  acoeff[7]= Ahat21;

  //  fprintf(stdout, "%f %f %f %f %f\n %f %f %f %f %f\n", Ares00, Ares10, Ares20, Ares30, Ares40, Ares01, Ares11, Ares21, Ares31, Ares41);
}

static void sf_cram_point3_fill_abins (sf_cram_point3 cram_point, float smp, 
                                       float ss, float sr, int ies, int ier,
                                       float oac, float ozc, float dac, float dzc,
                                       bool zoffset) {
    float dtw = DAM_TP*SF_PI/180.0,
          otw = OAM_TP*SF_PI/180.0; /* Tapering width in degrees */
    int dnw = 2*(int)(dtw/cram_point->db + 0.5), onw = (int)(otw/cram_point->db + 0.5); /* Taper width in samples */
    int isxy, jsxy, ihxy, jhxy, ia, iz, ida, ioa, ioz, idz;
    int ioaf, ioal, iozf, iozl, idaf, idal, idzf, idzl;
    float doa = 0.0, doz = 0.0, dda = 0.0, ddz = 0.0;
    float oa, da, oz, dz;
    float ds, dh, sa, sb, ha, hb, dw, tw = 1.0;
    bool flip = false;

    /* Integral weight */
    if (cram_point->amp) {
        /* See Bleistein et al (2003), equation (21) for details */
        if (cram_point->erefl) {
            dw = sqrtf (ss);
            dw /= sqrtf (fabsf (cram_point->src_exits[ies].j));
            dw *= cram_point->src_exits[ies].cs;
            sb = sinf (cram_point->b0 + cram_point->src_exits[ies].ib*cram_point->db);
            dw *= sqrtf (fabsf (sb));
            smp *= dw;
        } else {
            dw = sqrtf (ss*sr);
            dw *= cosf (0.5*oac);
            dw /= sqrtf (fabsf (cram_point->src_exits[ies].j*cram_point->rcv_exits[ier].j));
            dw *= cram_point->src_exits[ies].cs*cram_point->rcv_exits[ier].cs;
            sb = sinf (cram_point->b0 + cram_point->src_exits[ies].ib*cram_point->db);
            hb = zoffset ? sb : sinf (cram_point->b0 + cram_point->rcv_exits[ier].ib*cram_point->db);
            dw *= sqrtf (fabsf (sb*hb));
            smp *= dw;
        }
    }

    /* Compute tapeting weights at the edges of the mute zone */
    dw = cram_point->oam - fabsf (oac);
    if (dw < 2.0*otw) {
        if (dw < otw)
            tw = 0.0;
        else
            tw *= (dw - otw)/otw;
    }
    dw = cram_point->dam - fabsf (dac);
    if (dw < 2.0*dtw) {
        if (dw < dtw)
            tw = 0.0;
        else
            tw *= (dw - dtw)/dtw;
    }

    /* Contribute sample to the stacked image */
    cram_point->img += tw*smp;
    cram_point->hits += 1.0;

    if (!cram_point->agath && !cram_point->dipgath)
        return;

    if (cram_point->extrap && !zoffset) {
        /* Compute change in inclination and azimuth angles for the source and receiver branches
           using db/dx,db/dy and da/dx,da/dy with respect to a displacement away from the
           source and receiver on the surface; then find dip and scattering angle deviation,
           which correspond to this displacement */
        for (isxy = 0; isxy < 2; isxy++) { /* d/dx, d/dy, source side */
            for (jsxy = 0; jsxy < 2; jsxy++) { /* +/- shift in x/y on the source side */
                ds = jsxy != 0 ? cram_point->ds : -cram_point->ds;
                sb = cram_point->b0 +
                     (cram_point->src_exits[ies].ib + cram_point->src_exits[ies].ibxy[isxy]*ds)*
                     cram_point->db;
                sa = cram_point->a0 +
                     (cram_point->src_exits[ies].ia + cram_point->src_exits[ies].iaxy[isxy]*ds)*
                     cram_point->da;
                for (ihxy = 0; ihxy < 2; ihxy++) { /* d/dx, d/dy, receiver side */
                    for (jhxy = 0; jhxy < 2; jhxy++) { /* +/- shift in x/y on the receiver side */
                        dh = jhxy != 0 ? cram_point->ds : -cram_point->ds;
                        hb = cram_point->b0 +
                             (cram_point->rcv_exits[ier].ib + cram_point->rcv_exits[ier].ibxy[ihxy]*dh)*
                             cram_point->db;
                        ha = cram_point->a0 +
                             (cram_point->rcv_exits[ier].ia + cram_point->rcv_exits[ier].iaxy[ihxy]*dh)*
                             cram_point->da;
                        /* New system of subsurface angles for the surface displacment */
                        sf_cram_point3_angles (sb, sa, hb, ha, &oa, &oz, &da, &dz);
                        /* Scattering angle spread with respect to the initial values */
                        sf_cram_point3_aaz_spread (oac, ozc, oa, oz, &doa, &doz);
                        /* Dip angle spread with respect to the initial values */
                        sf_cram_point3_aaz_spread (dac, dzc, da, dz, &dda, &ddz);
                    }
                }
            }
        }
        /* Average changes in dip and scattering angles and their azimuths */
        doa /= 16.0;
        doz /= 16.0;
        dda /= 16.0;
        ddz /= 16.0;
        /* Extrapolate azimuth further away */
        doa *= 0.75;
        /* Extrapolate dip only in the vicinity of a receiver */
//.........这里部分代码省略.........

GLUSboolean update(GLUSfloat time)
{
    static GLfloat angle = 0.0f;

    GLuint primitivesWritten;

    // Field of view

    GLfloat rotationMatrix[16];

    GLfloat positionTextureSpace[4];
    GLfloat directionTextureSpace[3];
    GLfloat leftNormalTextureSpace[3];
    GLfloat rightNormalTextureSpace[3];
    GLfloat backNormalTextureSpace[3];

    GLfloat xzPosition2D[4];

    //

    GLfloat tmvpMatrix[16];

    // Animation update

    g_personView.cameraPosition[0] = -cosf(2.0f * PIf * angle / TURN_DURATION) * TURN_RADIUS * METERS_TO_VIRTUAL_WORLD_SCALE;
    g_personView.cameraPosition[2] = -sinf(2.0f * PIf * angle / TURN_DURATION) * TURN_RADIUS * METERS_TO_VIRTUAL_WORLD_SCALE;

    g_personView.cameraDirection[0] = sinf(2.0f * PIf * angle / TURN_DURATION);
    g_personView.cameraDirection[2] = -cosf(2.0f * PIf * angle / TURN_DURATION);

    if (g_animationOn)
    {
        angle += time;
    }

    glusLookAtf(g_viewMatrix, g_activeView->cameraPosition[0], g_activeView->cameraPosition[1], g_activeView->cameraPosition[2], g_activeView->cameraPosition[0] + g_activeView->cameraDirection[0], g_activeView->cameraPosition[1] + g_activeView->cameraDirection[1],
            g_activeView->cameraPosition[2] + g_activeView->cameraDirection[2], g_activeView->cameraUp[0], g_activeView->cameraUp[1], g_activeView->cameraUp[2]);

    glusMatrix4x4Identityf(tmvpMatrix);
    glusMatrix4x4Multiplyf(tmvpMatrix, tmvpMatrix, g_projectionMatrix);
    glusMatrix4x4Multiplyf(tmvpMatrix, tmvpMatrix, g_viewMatrix);
    glusMatrix4x4Multiplyf(tmvpMatrix, tmvpMatrix, g_textureToWorldMatrix);

    // Position

    xzPosition2D[0] = g_personView.cameraPosition[0];
    xzPosition2D[1] = 0.0f;
    xzPosition2D[2] = g_personView.cameraPosition[2];
    xzPosition2D[3] = g_personView.cameraPosition[3];

    glusMatrix4x4MultiplyPoint4f(positionTextureSpace, g_worldToTextureMatrix, xzPosition2D);

    // Direction

    glusMatrix4x4MultiplyVector3f(directionTextureSpace, g_worldToTextureMatrix, g_personView.cameraDirection);

    // Left normal of field of view

    glusMatrix4x4Identityf(rotationMatrix);
    glusMatrix4x4RotateRyf(rotationMatrix, g_personView.fov * (g_width / g_height) / 2.0f + 90.0f);
    glusMatrix4x4MultiplyVector3f(leftNormalTextureSpace, rotationMatrix, g_personView.cameraDirection);
    glusMatrix4x4MultiplyVector3f(leftNormalTextureSpace, g_worldToTextureNormalMatrix, leftNormalTextureSpace);

    // Right normal of field of view

    glusMatrix4x4Identityf(rotationMatrix);
    glusMatrix4x4RotateRyf(rotationMatrix, -g_personView.fov * (g_width / g_height) / 2.0f - 90.0f);
    glusMatrix4x4MultiplyVector3f(rightNormalTextureSpace, rotationMatrix, g_personView.cameraDirection);
    glusMatrix4x4MultiplyVector3f(rightNormalTextureSpace, g_worldToTextureNormalMatrix, rightNormalTextureSpace);

    // Back normal of field of view

    glusMatrix4x4Identityf(rotationMatrix);
    glusMatrix4x4RotateRyf(rotationMatrix, 180.0f);
    glusMatrix4x4MultiplyVector3f(backNormalTextureSpace, rotationMatrix, g_personView.cameraDirection);
    glusMatrix4x4MultiplyVector3f(backNormalTextureSpace, g_worldToTextureNormalMatrix, backNormalTextureSpace);

    // OpenGL stuff

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    // Pass one.

    // Disable any rasterization
    glEnable(GL_RASTERIZER_DISCARD);

    glUseProgram(g_programPassOne.program);

    glUniform4fv(g_positionTextureSpacePassOneLocation, 1, positionTextureSpace);
    glUniform3fv(g_leftNormalTextureSpacePassOneLocation, 1, leftNormalTextureSpace);
    glUniform3fv(g_rightNormalTextureSpacePassOneLocation, 1, rightNormalTextureSpace);
    glUniform3fv(g_backNormalTextureSpacePassOneLocation, 1, backNormalTextureSpace);

    glBindVertexArray(g_vaoPassOne);

    // Bind to vertices used in render pass two. To this buffer is written.
    glBindBufferBase(GL_TRANSFORM_FEEDBACK_BUFFER, 0, g_verticesPassTwoVBO);

    // We need to know, how many primitives are written. So start the query.
    glBeginQuery(GL_TRANSFORM_FEEDBACK_PRIMITIVES_WRITTEN, g_transformFeedbackQuery);
//.........这里部分代码省略.........

Array2D<float> evaluateImage(const Image &src, const std::vector<float> &scales, int steps)
{
    int border = 2.5f * (*std::max_element(scales.begin(), scales.end()));
    int w = src.getWidth(), h = src.getHeight();

    progress(0);

    Array2D<float> globalProfile[steps];
    for(int i=0; i<steps; i++) {
        globalProfile[i] = Array2D<float>(w,h);
        globalProfile[i].fill(1.f);
    }

    {
        Array2D<float> currentProfile[steps];
        for(int i=0; i<steps; i++)
            currentProfile[i] = Array2D<float>(w,h);
        
        PrefixSums ps(src);

        DifferenceJob jobs[steps];
        Completion c;

        float progress_step = 1.f/(steps * scales.size());
        float progress_done = 0.f;

        for(int i=0; i<steps; i++) {
            jobs[i].src = &ps;
            jobs[i].dst = &currentProfile[i];
            jobs[i].x1 = border;
            jobs[i].y1 = border;
            jobs[i].x2 = w-border;
            jobs[i].y2 = h-border;
            jobs[i].progress_step = progress_step;
            jobs[i].progress_done = &progress_done;
            jobs[i].completion = &c;
        }
        
        for(int s=0; s<(int)scales.size(); s++)
        {
            for(int i=0; i<steps; i++)
            {
                float angle = 2.0f*M_PI/steps*i;
                jobs[i].scale = scales[s];
                jobs[i].dx = jobs[i].scale * cosf(angle);
                jobs[i].dy = jobs[i].scale * sinf(angle);
                aq->queue(&jobs[i]);
            }
            
            c.wait();
            
            for(int i=0; i<steps; i++)
                for(int y=0; y<h; y++)
                    for(int x=0; x<w; x++)
                        globalProfile[i][y][x] *= currentProfile[i][y][x];
        }
    }

    Array2D<float> maxs(globalProfile[0]), mins(globalProfile[0]);
    
    for(int i=1; i<steps; i++)
        for(int y=border; y<h-border; y++)
            for(int x=border; x<w-border; x++) {
                maxs[y][x] = std::max(maxs[y][x], globalProfile[i][y][x]);
                mins[y][x] = std::min(mins[y][x], globalProfile[i][y][x]);
            }

    Array2D<float> eval(w,h);
    eval.fill(0.f);
    for(int y=border; y<h-border; y++)
        for(int x=border; x<w-border; x++)
            eval[y][x] = powf(maxs[y][x] - mins[y][x], 1.0f/scales.size());
    return eval;
}

void FlipX3D::update(float time)
{
    float angle = (float)M_PI * time; // 180 degrees
    float mz = sinf(angle);
    angle = angle / 2.0f; // x calculates degrees from 0 to 90
    float mx = cosf(angle);

    Vec3 v0, v1, v, diff;

    v0 = getOriginalVertex(Vec2(1.0f, 1.0f));
    v1 = getOriginalVertex(Vec2());

    float    x0 = v0.x;
    float    x1 = v1.x;
    float    x;
    Vec2    a, b, c, d;

    if ( x0 > x1 )
    {
        // Normal Grid
        a.setZero();
        b.set(0.0f, 1.0f);
        c.set(1.0f, 0.0f);
        d.set(1.0f, 1.0f);
        x = x0;
    }
    else
    {
        // Reversed Grid
        c.setZero();
        d.set(0.0f, 1.0f);
        a.set(1.0f, 0.0f);
        b.set(1.0f, 1.0f);
        x = x1;
    }
    
    diff.x = ( x - x * mx );
    diff.z = fabsf( floorf( (x * mz) / 4.0f ) );

    // bottom-left
    v = getOriginalVertex(a);
    v.x = diff.x;
    v.z += diff.z;
    setVertex(a, v);
    
    // upper-left
    v = getOriginalVertex(b);
    v.x = diff.x;
    v.z += diff.z;
    setVertex(b, v);
    
    // bottom-right
    v = getOriginalVertex(c);
    v.x -= diff.x;
    v.z -= diff.z;
    setVertex(c, v);
    
    // upper-right
    v = getOriginalVertex(d);
    v.x -= diff.x;
    v.z -= diff.z;
    setVertex(d, v);
}

/// update - update circle controller
void AC_Circle::update()
{
    // calculate dt
    uint32_t now = hal.scheduler->millis();
    float dt = (now - _last_update) / 1000.0f;

    // update circle position at 10hz
    if (dt > 0.095f) {

        // double check dt is reasonable
        if (dt >= 1.0f) {
            dt = 0.0;
        }
        // capture time since last iteration
        _last_update = now;

        // ramp up angular velocity to maximum
        if (_rate >= 0) {
            if (_angular_vel < _angular_vel_max) {
                _angular_vel += _angular_accel * dt;
                _angular_vel = constrain_float(_angular_vel, 0, _angular_vel_max);
            }
        } else {
            if (_angular_vel > _angular_vel_max) {
                _angular_vel += _angular_accel * dt;
                _angular_vel = constrain_float(_angular_vel, _angular_vel_max, 0);
            }
        }

        // update the target angle and total angle traveled
        float angle_change = _angular_vel * dt;
        _angle += angle_change;
        _angle = wrap_PI(_angle);
        _angle_total += angle_change;

        // if the circle_radius is zero we are doing panorama so no need to update loiter target
        if (_radius != 0.0) {
            // calculate target position
            Vector3f target;
            target.x = _center.x + _radius * cosf(-_angle);
            target.y = _center.y - _radius * sinf(-_angle);
            target.z = _pos_control.get_alt_target();

            // update position controller target
            _pos_control.set_pos_target(target);

            // heading is 180 deg from vehicles target position around circle
            _yaw = wrap_PI(_angle-PI) * AC_CIRCLE_DEGX100;
        } else {
            // set target position to center
            Vector3f target;
            target.x = _center.x;
            target.y = _center.y;
            target.z = _pos_control.get_alt_target();

            // update position controller target
            _pos_control.set_pos_target(target);

            // heading is same as _angle but converted to centi-degrees
            _yaw = _angle * AC_CIRCLE_DEGX100;
        }

        // trigger position controller on next update
        _pos_control.trigger_xy();
    }

    // run loiter's position to velocity step
    _pos_control.update_pos_controller(false);
}

//
/// \brief Generates geometry for a sphere.  Allocates memory for the vertex data and stores 
///        the results in the arrays.  Generate index list for a TRIANGLE_STRIP
/// \param numSlices The number of slices in the sphere
/// \param vertices If not NULL, will contain array of float3 positions
/// \param normals If not NULL, will contain array of float3 normals
/// \param texCoords If not NULL, will contain array of float2 texCoords
/// \param indices If not NULL, will contain the array of indices for the triangle strip
/// \return The number of indices required for rendering the buffers (the number of indices stored in the indices array
///         if it is not NULL ) as a GL_TRIANGLE_STRIP
//
int ESUTIL_API esGenSphere ( int numSlices, float radius, GLfloat **vertices, GLfloat **normals, 
                             GLfloat **texCoords, GLuint **indices )
{
   int i;
   int j;
   int numParallels = numSlices / 2;
   int numVertices = ( numParallels + 1 ) * ( numSlices + 1 );
   int numIndices = numParallels * numSlices * 6;
   float angleStep = (2.0f * ES_PI) / ((float) numSlices);

   // Allocate memory for buffers
   if ( vertices != NULL )
      *vertices = malloc ( sizeof(GLfloat) * 3 * numVertices );
   
   if ( normals != NULL )
      *normals = malloc ( sizeof(GLfloat) * 3 * numVertices );

   if ( texCoords != NULL )
      *texCoords = malloc ( sizeof(GLfloat) * 2 * numVertices );

   if ( indices != NULL )
      *indices = malloc ( sizeof(GLuint) * numIndices );

   for ( i = 0; i < numParallels + 1; i++ )
   {
      for ( j = 0; j < numSlices + 1; j++ )
      {
         int vertex = ( i * (numSlices + 1) + j ) * 3; 

         if ( vertices )
         {
            (*vertices)[vertex + 0] = radius * sinf ( angleStep * (float)i ) *
                                               sinf ( angleStep * (float)j );
            (*vertices)[vertex + 1] = radius * cosf ( angleStep * (float)i );
            (*vertices)[vertex + 2] = radius * sinf ( angleStep * (float)i ) *
                                               cosf ( angleStep * (float)j );
         }

         if ( normals )
         {
            (*normals)[vertex + 0] = (*vertices)[vertex + 0] / radius;
            (*normals)[vertex + 1] = (*vertices)[vertex + 1] / radius;
            (*normals)[vertex + 2] = (*vertices)[vertex + 2] / radius;
         }

         if ( texCoords )
         {
            int texIndex = ( i * (numSlices + 1) + j ) * 2;
            (*texCoords)[texIndex + 0] = (float) j / (float) numSlices;
            (*texCoords)[texIndex + 1] = ( 1.0f - (float) i ) / (float) (numParallels - 1 );
         }
      }
   }

   // Generate the indices
   if ( indices != NULL )
   {
      GLuint *indexBuf = (*indices);
      for ( i = 0; i < numParallels ; i++ ) 
      {
         for ( j = 0; j < numSlices; j++ )
         {
            *indexBuf++  = i * ( numSlices + 1 ) + j;
            *indexBuf++ = ( i + 1 ) * ( numSlices + 1 ) + j;
            *indexBuf++ = ( i + 1 ) * ( numSlices + 1 ) + ( j + 1 );

            *indexBuf++ = i * ( numSlices + 1 ) + j;
            *indexBuf++ = ( i + 1 ) * ( numSlices + 1 ) + ( j + 1 );
            *indexBuf++ = i * ( numSlices + 1 ) + ( j + 1 );
         }
      }
   }

   return numIndices;
}

void CShapeRenderer::DrawCone(const CMat4& mat, float fButtomRadius, float fTopRadius, float fHeight, CColor coneColor, CColor bottomColor, CColor topColor, bool bSolid) const
{
    DrawCircle(mat, fButtomRadius, bottomColor, bSolid);
    CMat4 translateMat;
    translateMat.SetTranslate(CVec3(0, fHeight, 0));
    CMat4 topMat = mat * translateMat;
    DrawCircle(topMat, fTopRadius, topColor, bSolid);
    if (bSolid)
    {
        CVec3 center(mat[12], mat[13], mat[14]);
        CVec3 upDirection = mat.GetUpVec3();
        upDirection = upDirection * fHeight;
        CVec3 topCenter = center + upDirection;
        CVertexPC point;
        point.position = topCenter;
        point.color = coneColor;
        std::vector<unsigned short> indicesData;
        std::vector<CVertexPC> vertexData;
        vertexData.push_back(point);
        static const float fStepRadians = DegreesToRadians(15);
        for (float fRadian = 0; fRadian <= MATH_PI_DOUBLE; fRadian += fStepRadians)
        {
            CVec3 pos(fButtomRadius * sinf(fRadian), 0, fButtomRadius * cosf(fRadian));
            pos *= mat;
            point.position = pos;
            vertexData.push_back(point);
            if (vertexData.size() >= 3)
            {
                unsigned short index = (unsigned short)vertexData.size();
                indicesData.push_back(0);
                indicesData.push_back(index - 2);
                indicesData.push_back(index - 1);
            }
        }
        CRenderManager::GetInstance()->RenderTriangle(vertexData, indicesData, true);
    }
    else
    {
        CVertexPC startPos, endPos;
        startPos.color = coneColor;
        endPos.color = coneColor;
        CVec3 topPoint(fTopRadius, fHeight, 0);
        CVec3 buttomPoint(fButtomRadius, 0, 0);
        startPos.position = topPoint * mat;
        endPos.position = buttomPoint * mat;
        CRenderManager::GetInstance()->RenderLine(startPos, endPos, 1.0f, true);
        topPoint.X() *= -1;
        buttomPoint.X() *= -1;
        startPos.position = topPoint * mat;
        endPos.position = buttomPoint * mat;
        CRenderManager::GetInstance()->RenderLine(startPos, endPos, 1.0f, true);
        topPoint = CVec3(0, fHeight, fTopRadius);
        buttomPoint = CVec3(0, 0, fButtomRadius);
        startPos.position = topPoint * mat;
        endPos.position = buttomPoint * mat;
        CRenderManager::GetInstance()->RenderLine(startPos, endPos, 1.0f, true);
        topPoint.Z() *= -1;
        buttomPoint.Z() *= -1;
        startPos.position = topPoint * mat;
        endPos.position = buttomPoint * mat;
        CRenderManager::GetInstance()->RenderLine(startPos, endPos, 1.0f, true);
    }
}

    static bool HandleWarpCommand(ChatHandler* handler, char const* args)
    {
        if (!*args)
            return false;

        Player* player = handler->GetSession()->GetPlayer();

        char* arg1 = strtok((char*)args, " ");
        char* arg2 = strtok(NULL, " ");

        if (!arg1 || !arg2)
        return false;

        char dir = arg1[0];
        float value = float(atof(arg2));
        float x = player->GetPositionX();
        float y = player->GetPositionY();
        float z = player->GetPositionZ();
        float o = player->GetOrientation();
		 uint32 mapid = player->GetMapId();
        Map const* map = sMapMgr->CreateBaseMap(mapid);
        z = std::max(map->GetHeight(x, y, MAX_HEIGHT), map->GetWaterLevel(x, y));

    switch (dir)
    {
	case 'l': // left
        {
			x = x + cos(o+(M_PI/2))*value;
           y = y + sin(o+(M_PI/2))*value;

           player->TeleportTo(mapid, x, y, z, o);
        }
        break;
    case 'r': // right
        {
			x = x + cos(o-(M_PI/2))*value;
           y = y + sin(o-(M_PI/2))*value;

           player->TeleportTo(mapid, x, y, z, o);
        }
        break;
    case 'f': // forward
        {
			x = x + cosf(o)*value;
           y = y + sinf(o)*value;

           player->TeleportTo(mapid, x, y, z, o);
        }
        break;
    case 'u': // up
		{
           player->TeleportTo(mapid, x, y, z + value, o);
		}
        break;
	 case 'd': // down
		{
           player->TeleportTo(mapid, x, y, z - value, o);
		}
        break;
    case 'o': //orientation
        {
			o = Position::NormalizeOrientation((value * M_PI_F/180.0f)+ o);

           player->TeleportTo(mapid, x, y, z, o);
        }
        break;
    }
    return true;
  };

/*
// Getting feature map for the selected subimage  
//
// API
// int getFeatureMaps(const IplImage * image, const int k, featureMap **map);
// INPUT
// image             - selected subimage
// k                 - size of cells
// OUTPUT
// map               - feature map
// RESULT
// Error status
*/
int getFeatureMaps_dp(const IplImage * image,const int k, CvLSVMFeatureMap **map)
{
    int sizeX, sizeY;
    int p, px, strsz;
    int height, width, channels;
    int i, j, kk, c, ii, jj, d;
    float  * datadx, * datady;
    float tmp, x, y, tx, ty;
    IplImage * dx, * dy;
    int *nearest_x, *nearest_y;
    float *w, a_x, b_x;

    float kernel[3] = {-1.f, 0.f, 1.f}; 
    CvMat kernel_dx = cvMat(1, 3, CV_32F, kernel);
    CvMat kernel_dy = cvMat(3, 1, CV_32F, kernel);

    float * r;
    int    * alfa;
    
    float boundary_x[CNTPARTION+1];
    float boundary_y[CNTPARTION+1];
    float max, tmp_scal;
    int    maxi;

	height = image->height;
	width  = image->width ;

    channels  = image->nChannels;

	dx    = cvCreateImage(cvSize(image->width , image->height) , IPL_DEPTH_32F , 3);
    dy    = cvCreateImage(cvSize(image->width , image->height) , IPL_DEPTH_32F , 3);

    sizeX = width  / k;
    sizeY = height / k;
    px    = CNTPARTION  + 2 * CNTPARTION; // контрастное и не контрастное изображение
    p     = px;
    strsz = sizeX * p;
    allocFeatureMapObject(map, sizeX, sizeY, p,  px);

	cvFilter2D(image, dx, &kernel_dx, cvPoint(-1, 0));
	cvFilter2D(image, dy, &kernel_dy, cvPoint(0, -1));
	
    for(i = 0; i <= CNTPARTION; i++)
    {
        boundary_x[i] = cosf((((float)i) * (((float)PI) / (float) (CNTPARTION))));
        boundary_y[i] = sinf((((float)i) * (((float)PI) / (float) (CNTPARTION))));
    }/*for(i = 0; i <= CNTPARTION; i++) */

    r    = (float *)malloc( sizeof(float) * (width * height));
    alfa = (int   *)malloc( sizeof(int  ) * (width * height * 2));

    for(j = 1; j < height-1; j++)
    {
        datadx = (float*)(dx->imageData + dx->widthStep *j);
        datady = (float*)(dy->imageData + dy->widthStep *j);
        for(i = 1; i < width-1; i++)
        {
			c = 0;
            x = (datadx[i*channels+c]);
            y = (datady[i*channels+c]);

            r[j * width + i] =sqrtf(x*x + y*y);
            for(kk = 1; kk < channels; kk++)
            {
                tx = (datadx[i*channels+kk]);
                ty = (datady[i*channels+kk]);
                tmp =sqrtf(tx*tx + ty*ty);
                if(tmp > r[j * width + i])
                {
                    r[j * width + i] = tmp;
                    c = kk;
                    x = tx;
                    y = ty;
                }
            }/*for(kk = 1; kk < channels; kk++)*/

            
            
            max  = boundary_x[0]*x + boundary_y[0]*y;
            maxi = 0;
            for (kk = 0; kk < CNTPARTION; kk++) {
                tmp_scal = boundary_x[kk]*x + boundary_y[kk]*y;
                if (tmp_scal> max) {
                    max = tmp_scal;
                    maxi = kk;
                }else if (-tmp_scal> max) {
                    max = -tmp_scal;
//.........这里部分代码省略.........

// Determine if learning of wind and magnetic field will be enabled and set corresponding indexing limits to
// avoid unnecessary operations
void NavEKF2_core::setWindMagStateLearningMode()
{
    // If we are on ground, or in constant position mode, or don't have the right vehicle and sensing to estimate wind, inhibit wind states
    bool setWindInhibit = (!useAirspeed() && !assume_zero_sideslip()) || onGround || (PV_AidingMode == AID_NONE);
    if (!inhibitWindStates && setWindInhibit) {
        inhibitWindStates = true;
    } else if (inhibitWindStates && !setWindInhibit) {
        inhibitWindStates = false;
        // set states and variances
        if (yawAlignComplete && useAirspeed()) {
            // if we have airspeed and a valid heading, set the wind states to the reciprocal of the vehicle heading
            // which assumes the vehicle has launched into the wind
             Vector3f tempEuler;
            stateStruct.quat.to_euler(tempEuler.x, tempEuler.y, tempEuler.z);
            float windSpeed =  sqrtf(sq(stateStruct.velocity.x) + sq(stateStruct.velocity.y)) - tasDataDelayed.tas;
            stateStruct.wind_vel.x = windSpeed * cosf(tempEuler.z);
            stateStruct.wind_vel.y = windSpeed * sinf(tempEuler.z);

            // set the wind sate variances to the measurement uncertainty
            for (uint8_t index=22; index<=23; index++) {
                P[index][index] = sq(constrain_float(frontend->_easNoise, 0.5f, 5.0f) * constrain_float(_ahrs->get_EAS2TAS(), 0.9f, 10.0f));
            }
        } else {
            // set the variances using a typical wind speed
            for (uint8_t index=22; index<=23; index++) {
                P[index][index] = sq(5.0f);
            }
        }
    }

    // determine if the vehicle is manoevring
    if (accNavMagHoriz > 0.5f) {
        manoeuvring = true;
    } else {
        manoeuvring = false;
    }

    // Determine if learning of magnetic field states has been requested by the user
    uint8_t magCal = effective_magCal();
    bool magCalRequested =
            ((magCal == 0) && inFlight) || // when flying
            ((magCal == 1) && manoeuvring)  || // when manoeuvring
            ((magCal == 3) && firstInflightYawInit && firstInflightMagInit) || // when initial in-air yaw and mag field reset is complete
            (magCal == 4); // all the time

    // Deny mag calibration request if we aren't using the compass, it has been inhibited by the user,
    // we do not have an absolute position reference or are on the ground (unless explicitly requested by the user)
    bool magCalDenied = !use_compass() || (magCal == 2) || (onGround && magCal != 4);

    // Inhibit the magnetic field calibration if not requested or denied
    bool setMagInhibit = !magCalRequested || magCalDenied;
    if (!inhibitMagStates && setMagInhibit) {
        inhibitMagStates = true;
    } else if (inhibitMagStates && !setMagInhibit) {
        inhibitMagStates = false;
        // when commencing use of magnetic field states, set the variances equal to the observation uncertainty
        for (uint8_t index=16; index<=21; index++) {
            P[index][index] = sq(frontend->_magNoise);
        }
        // request a reset of the yaw and magnetic field states if not done before
        if (!magStateInitComplete || (!firstInflightMagInit && inFlight)) {
            magYawResetRequest = true;
        }
    }

    // If on ground we clear the flag indicating that the magnetic field in-flight initialisation has been completed
    // because we want it re-done for each takeoff
    if (onGround) {
        firstInflightYawInit = false;
        firstInflightMagInit = false;
    }

    // Adjust the indexing limits used to address the covariance, states and other EKF arrays to avoid unnecessary operations
    // if we are not using those states
    if (inhibitMagStates && inhibitWindStates) {
        stateIndexLim = 15;
    } else if (inhibitWindStates) {
        stateIndexLim = 21;
    } else {
        stateIndexLim = 23;
    }
}

/**
 * body()
 *
 * Creates a display list for the body of a person.
 */
void body()
{
   GLfloat fan_angle;   //angle for a triangle fan
   int sections = 40;   //number of triangles in a triangle fan

   //generate a display list for a body
   body_list_id = glGenLists(1);
   //compile the new display list
   glNewList(body_list_id, GL_COMPILE);

   glPushMatrix();
      glTranslatef(0.0f, 5.05f, 0.0f);
      glRotatef(90, 1.0, 0.0, 0.0);

      //create a torso with a cylinder
      glPushMatrix();
         GLUquadricObj * torso;
         torso = gluNewQuadric();
         gluQuadricDrawStyle(torso, GLU_FILL);
         gluCylinder(torso, 1.0f, 1.0f, 2.25f, 25, 25);
      glPopMatrix();

      //triangle fan for top of body
      glPushMatrix();
         glRotatef(180.0, 0.0, 1.0, 0.0);
         glBegin(GL_TRIANGLE_FAN);
         glNormal3f(0.0f, 0.0f, 1.0f);
         glVertex3f(0.0f, 0.0f, 0.0f);
         for(int i = 0; i <= sections; ++i) {
            fan_angle = (GLfloat) i * 2 * PI / sections;
            glVertex3f(cosf(fan_angle), sinf(fan_angle), 0.0f);
         }
         glEnd();
      glPopMatrix();

      //triangle fan for bottom of torso
      glPushMatrix();
         glTranslatef(0.0, 0.0, 2.25);
         glBegin(GL_TRIANGLE_FAN);
         glNormal3f(0.0f, 0.0f, 1.0f);
         glVertex3f(0.0f, 0.0f, 0.0f);
         for(int i = 0; i <= sections; ++i) {
            fan_angle = (GLfloat) i * 2 * PI / sections;
            glVertex3f(cosf(fan_angle), sinf(fan_angle), 0.0f);
         }
         glEnd();
      glPopMatrix();

      //arm 1
      glPushMatrix();
         glTranslatef(1.25f, 0.0f, 0.0f);
         GLUquadricObj * arm1;
         arm1 = gluNewQuadric();
         gluQuadricDrawStyle(arm1, GLU_FILL);
         gluCylinder(arm1, 0.25f, 0.25f, 2.5f, 12, 12);
      glPopMatrix();

      //triangle fan for top of arm 1
      glPushMatrix();
         glTranslatef(1.25, 0.0, 0.0);
         glRotatef(180.0, 0.0, 1.0, 0.0);
         glBegin(GL_TRIANGLE_FAN);
         glNormal3f(0.0f, 0.0f, 1.0f);
         glVertex3f(0.0f, 0.0f, 0.0f);
         sections = 12;
         for(int i = 0; i <= sections; ++i) {
            fan_angle = (GLfloat) i * 2 * PI / sections;
            glVertex3f(0.25*cosf(fan_angle), 0.25*sinf(fan_angle), 0.0f);
         }
         glEnd();
      glPopMatrix();

      //triangle fan for bottom of arm 1
      glPushMatrix();
         glTranslatef(1.25, 0.0, 2.5);
         glBegin(GL_TRIANGLE_FAN);
         glNormal3f(0.0f, 0.0f, 1.0f);
         glVertex3f(0.0f, 0.0f, 0.0f);
         sections = 12;
         for(int i = 0; i <= sections; ++i) {
            fan_angle = (GLfloat) i * 2 * PI / sections;
            glVertex3f(0.25*cosf(fan_angle), 0.25*sinf(fan_angle), 0.0f);
         }
         glEnd();
      glPopMatrix();

      //arm 2
      glPushMatrix();
         glTranslatef(-1.25f, 0.0f, 0.0f);
         GLUquadricObj * arm2;
         arm2 = gluNewQuadric();
         gluQuadricDrawStyle(arm2, GLU_FILL);
         gluCylinder(arm2, 0.25f, 0.25f, 2.5f, 12, 12);
      glPopMatrix();

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

void draw(HyperspaceSaverSettings *inSettings){
	if(inSettings->first){
		if(inSettings->dUseTunnels){  // only tunnels use caustic textures
			glDisable(GL_FOG);
			// Caustic textures can only be created after rendering context has been created
			// because they have to be drawn and then read back from the framebuffer.
#ifdef WIN32
			if(doingPreview) // super fast for Windows previewer
				inSettings->theCausticTextures = new causticTextures(8, inSettings->numAnimTexFrames, 32, 32, 1.0f, 0.01f, 10.0f);
			else  // normal
#endif
				inSettings->theCausticTextures = new causticTextures(8, inSettings->numAnimTexFrames, 100, 256, 1.0f, 0.01f, 20.0f);
			glEnable(GL_FOG);
		}
		if(inSettings->dShaders){
#ifdef WIN32
			if(doingPreview) // super fast for Windows previewer
				inSettings->theWNCM = new wavyNormalCubeMaps(inSettings->numAnimTexFrames, 32);
			else  // normal
#endif
				inSettings->theWNCM = new wavyNormalCubeMaps(inSettings->numAnimTexFrames, 128);
		}
		glViewport(inSettings->viewport[0], inSettings->viewport[1], inSettings->viewport[2], inSettings->viewport[3]);
		inSettings->first = 0;
	}

	// Variables for printing text
	static float computeTime = 0.0f;
	static float drawTime = 0.0f;
	//static rsTimer computeTimer, drawTimer;
	// start compute time timer
	//computeTimer.tick();

	glMatrixMode(GL_MODELVIEW);

	// Camera movements
	static float camHeading[3] = {0.0f, 0.0f, 0.0f};  // current, target, and last
	static int changeCamHeading = 1;
	static float camHeadingChangeTime[2] = {20.0f, 0.0f};  // total, elapsed
	static float camRoll[3] = {0.0f, 0.0f, 0.0f};  // current, target, and last
	static int changeCamRoll = 1;
	static float camRollChangeTime[2] = {1.0f, 0.0f};  // total, elapsed
	camHeadingChangeTime[1] += inSettings->frameTime;
	if(camHeadingChangeTime[1] >= camHeadingChangeTime[0]){  // Choose new direction
		camHeadingChangeTime[0] = rsRandf(15.0f) + 5.0f;
		camHeadingChangeTime[1] = 0.0f;
		camHeading[2] = camHeading[1];  // last = target
		if(changeCamHeading){
			// face forward most of the time
			if(rsRandi(6))
				camHeading[1] = 0.0f;
			// face backward the rest of the time
			else if(rsRandi(2))
				camHeading[1] = RS_PI;
			else
				camHeading[1] = -RS_PI;
			changeCamHeading = 0;
		}
		else
			changeCamHeading = 1;
	}
	float t = camHeadingChangeTime[1] / camHeadingChangeTime[0];
	t = 0.5f * (1.0f - cosf(RS_PI * t));
	camHeading[0] = camHeading[1] * t + camHeading[2] * (1.0f - t);
	camRollChangeTime[1] += inSettings->frameTime;
	if(camRollChangeTime[1] >= camRollChangeTime[0]){  // Choose new roll angle
		camRollChangeTime[0] = rsRandf(5.0f) + 10.0f;
		camRollChangeTime[1] = 0.0f;
		camRoll[2] = camRoll[1];  // last = target
		if(changeCamRoll){
			camRoll[1] = rsRandf(RS_PIx2*2) - RS_PIx2;
			changeCamRoll = 0;
		}
		else
			changeCamRoll = 1;
	}
	t = camRollChangeTime[1] / camRollChangeTime[0];
	t = 0.5f * (1.0f - cosf(RS_PI * t));
	camRoll[0] = camRoll[1] * t + camRoll[2] * (1.0f - t);

	static float pathDir[3] = {0.0f, 0.0f, -1.0f};
	inSettings->thePath->moveAlongPath(float(inSettings->dSpeed) * inSettings->frameTime * 0.04f);
	inSettings->thePath->update(inSettings->frameTime);
	inSettings->thePath->getPoint(inSettings->dDepth + 2, inSettings->thePath->step, inSettings->camPos);
	inSettings->thePath->getBaseDirection(inSettings->dDepth + 2, inSettings->thePath->step, pathDir);
	float pathAngle = atan2f(-pathDir[0], -pathDir[2]);

	glLoadIdentity();
	glRotatef((pathAngle + camHeading[0]) * RS_RAD2DEG, 0, 1, 0);
	glRotatef(camRoll[0] * RS_RAD2DEG, 0, 0, 1);
	glGetFloatv(GL_MODELVIEW_MATRIX, inSettings->billboardMat);
	glLoadIdentity();
	glRotatef(-camRoll[0] * RS_RAD2DEG, 0, 0, 1);
	glRotatef((-pathAngle - camHeading[0]) * RS_RAD2DEG, 0, 1, 0);
	glTranslatef(inSettings->camPos[0]*-1, inSettings->camPos[1]*-1, inSettings->camPos[2]*-1);
	glGetDoublev(GL_MODELVIEW_MATRIX, inSettings->modelMat);
	inSettings->unroll = camRoll[0] * RS_RAD2DEG;

	if(inSettings->dUseGoo){
		// calculate diagonal fov
//.........这里部分代码省略.........

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

    /*JUMP命令 */
    else if (iname == 39) {
      if (ist.name[1] == -100) {
	freg[2] = atanf(freg[2]);
	nextpc = pop(&call_stack) - 1;
      } else if (ist.name[1] == -200) {
	freg[2] = sqrtf(freg[2]);
	nextpc = pop(&call_stack) - 1;
      } else if (ist.name[1] == -300) {
	regist[4] = read_int();
	printf("read_int: %d\n", regist[4]);
	nextpc = pop(&call_stack) - 1;
      } else if (ist.name[1] == -400) {
	freg[2] = read_float();
	printf("read_float: %f\n", freg[2]);
	nextpc = pop(&call_stack) - 1;
      } else if (ist.name[1] == -500) {
	printf("aaa\n");
	if(ppp < freg[2])
	  {
	    ppp = freg[2];
	  }
	if(nnn > freg[2])
	  {
	    nnn = freg[2];
	  }
	//print_register();
	nextpc = label_trans(program->insts[pc].name[1], program) - 1;
	//freg[2] = sinf(freg[2]);
	//printf("read_float: %f\n", freg[2]);
	//nextpc = pop(&call_stack) - 1;
      } else if (ist.name[1] == -600) {
	freg[2] = cosf(freg[2]);
	//printf("read_float: %f\n", freg[2]);
	nextpc = pop(&call_stack) - 1;
      } else if (ist.name[1] == -700) {
	//print_int
	print_int();
	nextpc = pop(&call_stack) - 1;
      } else if (ist.name[1] == -800) {
	print_float();
	nextpc = pop(&call_stack) - 1;
      }
      else
	{
	  nextpc = ist.name[1] - 1;
	}
      /*
	if(strncmp(program->insts[pc].name[1] , "L_main", 6) == 0)
	{
	print_memory();
	}
      */
    }


    else if (iname == 40) {
          
      nextpc = ist.name[1] - 1;
      /*特殊ラベル*/
      if (ist.name[1] == -100) {
	freg[2] = atanf(freg[2]);
	nextpc = pc;
      } else if (ist.name[1] == -200) {
	freg[2] = sqrtf(freg[2]);

int main( int argc, char* argv[] )
{
    // Length of vectors
    unsigned int n = 100000;
 
    // Host input vectors
    double *h_a;
    double *h_b;
    // Host output vector
    double *h_c;
 
    // Device input buffers
    cl_mem d_a;
    cl_mem d_b;
    // Device output buffer
    cl_mem d_c;
 
    cl_platform_id cpPlatform;        // OpenCL platform
    cl_device_id device_id;           // device ID
    cl_context context;               // context
    cl_command_queue queue;           // command queue
    cl_program program;               // program
    cl_kernel kernel;                 // kernel
 
    // Size, in bytes, of each vector
    size_t bytes = n*sizeof(double);
 
    // Allocate memory for each vector on host
    h_a = (double*)malloc(bytes);
    h_b = (double*)malloc(bytes);
    h_c = (double*)malloc(bytes);
 
    // Initialize vectors on host
    int i;
    for( i = 0; i < n; i++ )
    {
        h_a[i] = sinf(i)*sinf(i);
        h_b[i] = cosf(i)*cosf(i);
    }
 
    size_t globalSize, localSize;
    cl_int err;
 
    // Number of work items in each local work group
    localSize = 64;
 
    // Number of total work items - localSize must be devisor
    globalSize = ceil(n/(float)localSize)*localSize;
 
    // Bind to platform
    err = clGetPlatformIDs(1, &cpPlatform, NULL);
 
    // Get ID for the device
    err = clGetDeviceIDs(cpPlatform, CL_DEVICE_TYPE_GPU, 1, &device_id, NULL);
 
    // Create a context  
    context = clCreateContext(0, 1, &device_id, NULL, NULL, &err);
 
    // Create a command queue 
    queue = clCreateCommandQueue(context, device_id, 0, &err);
 
    // Create the compute program from the source buffer
    program = clCreateProgramWithSource(context, 1,
                            (const char **) & kernelSource, NULL, &err);
 
    // Build the program executable 
    clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
 
    // Create the compute kernel in the program we wish to run
    kernel = clCreateKernel(program, "vecAdd", &err);
 
    // Create the input and output arrays in device memory for our calculation
    d_a = clCreateBuffer(context, CL_MEM_READ_ONLY, bytes, NULL, NULL);
    d_b = clCreateBuffer(context, CL_MEM_READ_ONLY, bytes, NULL, NULL);
    d_c = clCreateBuffer(context, CL_MEM_WRITE_ONLY, bytes, NULL, NULL);
 
    // Write our data set into the input array in device memory
    err = clEnqueueWriteBuffer(queue, d_a, CL_TRUE, 0,
                                   bytes, h_a, 0, NULL, NULL);
    err |= clEnqueueWriteBuffer(queue, d_b, CL_TRUE, 0,
                                   bytes, h_b, 0, NULL, NULL);
 
    // Set the arguments to our compute kernel
    err  = clSetKernelArg(kernel, 0, sizeof(cl_mem), &d_a);
    err |= clSetKernelArg(kernel, 1, sizeof(cl_mem), &d_b);
    err |= clSetKernelArg(kernel, 2, sizeof(cl_mem), &d_c);
    err |= clSetKernelArg(kernel, 3, sizeof(unsigned int), &n);
 
    // Execute the kernel over the entire range of the data set  
    err = clEnqueueNDRangeKernel(queue, kernel, 1, NULL, &globalSize, &localSize,
                                                              0, NULL, NULL);
 
    // Wait for the command queue to get serviced before reading back results
    clFinish(queue);
 
    // Read the results from the device
    clEnqueueReadBuffer(queue, d_c, CL_TRUE, 0,
                                bytes, h_c, 0, NULL, NULL );
 
    //Sum up vector c and print result divided by n, this should equal 1 within error
//.........这里部分代码省略.........