//.........这里部分代码省略.........
                                    &paohSDERasterBands, 
                                    (long*)&nBands);
    if( nSDEErr != SE_SUCCESS )
    {
        IssueSDEError( nSDEErr, "SE_raster_get_bands" );
        return CE_Fatal;
    }
    
    SE_RASBANDINFO band;
    
    // grab our other stuff from the first band and hope for the best
    band = paohSDERasterBands[0];
    
    
    nSDEErr = SE_rasbandinfo_get_band_size( band, 
                                            (long*)&nRasterXSize, 
                                            (long*)&nRasterYSize);
    if( nSDEErr != SE_SUCCESS )
    {
        IssueSDEError( nSDEErr, "SE_rasbandinfo_get_band_size" );
        return CE_Fatal;
    }
    
    SE_ENVELOPE extent;
    nSDEErr = SE_rasbandinfo_get_extent(band, &extent);
    if( nSDEErr != SE_SUCCESS )
    {
        IssueSDEError( nSDEErr, "SE_rasbandinfo_get_extent" );
        return CE_Fatal;
    }
    dfMinX = extent.minx;
    dfMinY = extent.miny;
    dfMaxX = extent.maxx;
    dfMaxY = extent.maxy;

    CPLDebug("SDERASTER", "minx: %.5f, miny: %.5f, maxx: %.5f, maxy: %.5f", dfMinX, dfMinY, dfMaxX, dfMaxY);

    // x0 roughly corresponds to dfMinX
    // y0 roughly corresponds to dfMaxY
    // They can be different than the extent parameters because 
    // SDE uses offsets.  The following info is from Duarte Carreira 
    // in relation to bug #2063 http://trac.osgeo.org/gdal/ticket/2063 :

    // Depending on how the data was loaded, the pixel width 
    // or pixel height may include a pixel offset from the nearest
    // tile boundary. An offset will be indicated by aplus sign 
    // "+" followed by a value. The value indicates the number 
    // of pixels the nearest tile boundary is to the left of
    // the image for the x dimension or the number of
    // pixels the nearest tile boundary is above the image for 
    // the y dimension. The offset information is only useful 
    // for advanced application developers who need to know 
    // where the image begins in relation to the underlying tile structure

    LFLOAT x0, y0; 
    nSDEErr = SE_rasbandinfo_get_tile_origin(band, &x0, &y0);
    if( nSDEErr != SE_SUCCESS )
    {
        IssueSDEError( nSDEErr, "SE_rasbandinfo_get_tile_origin" );
        return CE_Fatal;
    }
    CPLDebug("SDERASTER", "Tile origin: %.5f, %.5f", x0, y0);

    // we also need to adjust dfMaxX and dfMinY otherwise the cell size will change
    dfMaxX = (x0-dfMinX) + dfMaxX;
    dfMinY = (y0-dfMaxY) + dfMinY;

    // adjust the bbox based on the tile origin.
    dfMinX = MIN(x0, dfMinX);
    dfMaxY = MAX(y0, dfMaxY);

    nSDEErr = SE_rasterattr_create(&hAttributes, false);
    if( nSDEErr != SE_SUCCESS )
    {
        IssueSDEError( nSDEErr, "SE_rasterattr_create" );
        return CE_Fatal;
    }
    
    // Grab the pointer for our member variable

    nSDEErr = SE_stream_create(hConnection, &hStream);
    if( nSDEErr != SE_SUCCESS )
    {
        IssueSDEError( nSDEErr, "SE_stream_create" );
        return CE_Fatal;
    }

    
    for (int i=0; i < nBands; i++) {
        SetBand( i+1, new SDERasterBand( this, i+1, -1, &(paohSDERasterBands[i]) ));
    }

    GDALRasterBand* b = GetRasterBand(1);
    
    eDataType = b->GetRasterDataType();
    
    SE_rasterinfo_free(raster);

    return CE_None;
}

initial_setup (j_decompress_ptr cinfo)
/* Called once, when first SOS marker is reached */
{
    int ci;
    jpeg_component_info *compptr;

    /* Make sure image isn't bigger than I can handle */
    if ((long) cinfo->image_height > (long) JPEG_MAX_DIMENSION ||
            (long) cinfo->image_width > (long) JPEG_MAX_DIMENSION)
        ERREXIT1(cinfo, JERR_IMAGE_TOO_BIG, (unsigned int) JPEG_MAX_DIMENSION);

    if (cinfo->process == JPROC_LOSSLESS) {
        /* If precision > compiled-in value, we must downscale */
        if (cinfo->data_precision > BITS_IN_JSAMPLE)
            WARNMS2(cinfo, JWRN_MUST_DOWNSCALE,
                    cinfo->data_precision, BITS_IN_JSAMPLE);
    }
    else {  /* Lossy processes */
        /* For now, precision must match compiled-in value... */
        if (cinfo->data_precision != BITS_IN_JSAMPLE)
            ERREXIT1(cinfo, JERR_BAD_PRECISION, cinfo->data_precision);
    }

    /* Check that number of components won't exceed internal array sizes */
    if (cinfo->num_components > MAX_COMPONENTS)
        ERREXIT2(cinfo, JERR_COMPONENT_COUNT, cinfo->num_components,
                 MAX_COMPONENTS);

    /* Compute maximum sampling factors; check factor validity */
    cinfo->max_h_samp_factor = 1;
    cinfo->max_v_samp_factor = 1;
    for (ci = 0, compptr = cinfo->comp_info; ci < cinfo->num_components;
            ci++, compptr++) {
        if (compptr->h_samp_factor<=0 || compptr->h_samp_factor>MAX_SAMP_FACTOR ||
                compptr->v_samp_factor<=0 || compptr->v_samp_factor>MAX_SAMP_FACTOR)
            ERREXIT(cinfo, JERR_BAD_SAMPLING);
        cinfo->max_h_samp_factor = MAX(cinfo->max_h_samp_factor,
                                       compptr->h_samp_factor);
        cinfo->max_v_samp_factor = MAX(cinfo->max_v_samp_factor,
                                       compptr->v_samp_factor);
    }

    /* We initialize codec_data_unit and min_codec_data_unit to data_unit.
     * In the full decompressor, this will be overridden by jdmaster.c;
     * but in the transcoder, jdmaster.c is not used, so we must do it here.
     */
    cinfo->min_codec_data_unit = cinfo->data_unit;

    /* Compute dimensions of components */
    for (ci = 0, compptr = cinfo->comp_info; ci < cinfo->num_components;
            ci++, compptr++) {
        compptr->codec_data_unit = cinfo->data_unit;
        /* Size in data units */
        compptr->width_in_data_units = (JDIMENSION)
                                       jdiv_round_up((long) cinfo->image_width * (long) compptr->h_samp_factor,
                                               (long) (cinfo->max_h_samp_factor * cinfo->data_unit));
        compptr->height_in_data_units = (JDIMENSION)
                                        jdiv_round_up((long) cinfo->image_height * (long) compptr->v_samp_factor,
                                                (long) (cinfo->max_v_samp_factor * cinfo->data_unit));
        /* downsampled_width and downsampled_height will also be overridden by
         * jdmaster.c if we are doing full decompression.  The transcoder library
         * doesn't use these values, but the calling application might.
         */
        /* Size in samples */
        compptr->downsampled_width = (JDIMENSION)
                                     jdiv_round_up((long) cinfo->image_width * (long) compptr->h_samp_factor,
                                             (long) cinfo->max_h_samp_factor);
        compptr->downsampled_height = (JDIMENSION)
                                      jdiv_round_up((long) cinfo->image_height * (long) compptr->v_samp_factor,
                                              (long) cinfo->max_v_samp_factor);
        /* Mark component needed, until color conversion says otherwise */
        compptr->component_needed = TRUE;
        /* Mark no quantization table yet saved for component */
        compptr->quant_table = NULL;
    }

    /* Compute number of fully interleaved MCU rows. */
    cinfo->total_iMCU_rows = (JDIMENSION)
                             jdiv_round_up((long) cinfo->image_height,
                                           (long) (cinfo->max_v_samp_factor*cinfo->data_unit));

    /* Decide whether file contains multiple scans */
    if (cinfo->comps_in_scan < cinfo->num_components ||
            cinfo->process == JPROC_PROGRESSIVE)
        cinfo->inputctl->has_multiple_scans = TRUE;
    else
        cinfo->inputctl->has_multiple_scans = FALSE;
}

// host stub function
void ops_par_loop_update_halo_kernel2_zvel_plus_4_left(char const *name, ops_block block, int dim, int* range,
 ops_arg arg0, ops_arg arg1, ops_arg arg2) {
  ops_arg args[3] = { arg0, arg1, arg2};


  ops_timing_realloc(81,"update_halo_kernel2_zvel_plus_4_left");
  OPS_kernels[81].count++;

  //compute locally allocated range for the sub-block
  int start[3];
  int end[3];
  #ifdef OPS_MPI
  sub_block_list sb = OPS_sub_block_list[block->index];
  if (!sb->owned) return;
  for ( int n=0; n<3; n++ ){
    start[n] = sb->decomp_disp[n];end[n] = sb->decomp_disp[n]+sb->decomp_size[n];
    if (start[n] >= range[2*n]) {
      start[n] = 0;
    }
    else {
      start[n] = range[2*n] - start[n];
    }
    if (sb->id_m[n]==MPI_PROC_NULL && range[2*n] < 0) start[n] = range[2*n];
    if (end[n] >= range[2*n+1]) {
      end[n] = range[2*n+1] - sb->decomp_disp[n];
    }
    else {
      end[n] = sb->decomp_size[n];
    }
    if (sb->id_p[n]==MPI_PROC_NULL && (range[2*n+1] > sb->decomp_disp[n]+sb->decomp_size[n]))
      end[n] += (range[2*n+1]-sb->decomp_disp[n]-sb->decomp_size[n]);
  }
  #else //OPS_MPI
  for ( int n=0; n<3; n++ ){
    start[n] = range[2*n];end[n] = range[2*n+1];
  }
  #endif //OPS_MPI

  int x_size = MAX(0,end[0]-start[0]);
  int y_size = MAX(0,end[1]-start[1]);
  int z_size = MAX(0,end[2]-start[2]);


  int xdim0 = args[0].dat->size[0]*args[0].dat->dim;
  int ydim0 = args[0].dat->size[1];
  int xdim1 = args[1].dat->size[0]*args[1].dat->dim;
  int ydim1 = args[1].dat->size[1];

  //build opencl kernel if not already built

  buildOpenCLKernels_update_halo_kernel2_zvel_plus_4_left(
  xdim0,ydim0,xdim1,ydim1);

  //Timing
  double t1,t2,c1,c2;
  ops_timers_core(&c2,&t2);

  //set up OpenCL thread blocks
  size_t globalWorkSize[3] = {((x_size-1)/OPS_block_size_x+ 1)*OPS_block_size_x, ((y_size-1)/OPS_block_size_y + 1)*OPS_block_size_y, MAX(1,end[2]-start[2])};
  size_t localWorkSize[3] =  {OPS_block_size_x,OPS_block_size_y,1};


  int *arg2h = (int *)arg2.data;

  int consts_bytes = 0;

  consts_bytes += ROUND_UP(NUM_FIELDS*sizeof(int));

  reallocConstArrays(consts_bytes);

  consts_bytes = 0;
  arg2.data = OPS_consts_h + consts_bytes;
  arg2.data_d = OPS_consts_d + consts_bytes;
  for (int d=0; d<NUM_FIELDS; d++) ((int *)arg2.data)[d] = arg2h[d];
  consts_bytes += ROUND_UP(NUM_FIELDS*sizeof(int));
  mvConstArraysToDevice(consts_bytes);
  int dat0 = args[0].dat->elem_size;
  int dat1 = args[1].dat->elem_size;

  //set up initial pointers
  int d_m[OPS_MAX_DIM];
  #ifdef OPS_MPI
  for (int d = 0; d < dim; d++) d_m[d] = args[0].dat->d_m[d] + OPS_sub_dat_list[args[0].dat->index]->d_im[d];
  #else //OPS_MPI
  for (int d = 0; d < dim; d++) d_m[d] = args[0].dat->d_m[d];
  #endif //OPS_MPI
  int base0 = 1 * 
  (start[0] * args[0].stencil->stride[0] - args[0].dat->base[0] - d_m[0]);
  base0 = base0 + args[0].dat->size[0] *
  (start[1] * args[0].stencil->stride[1] - args[0].dat->base[1] - d_m[1]);
  base0 = base0 + args[0].dat->size[0] *  args[0].dat->size[1] *
  (start[2] * args[0].stencil->stride[2] - args[0].dat->base[2] - d_m[2]);

  #ifdef OPS_MPI
  for (int d = 0; d < dim; d++) d_m[d] = args[1].dat->d_m[d] + OPS_sub_dat_list[args[1].dat->index]->d_im[d];
  #else //OPS_MPI
  for (int d = 0; d < dim; d++) d_m[d] = args[1].dat->d_m[d];
  #endif //OPS_MPI
  int base1 = 1 * 
//.........这里部分代码省略.........

lapack_int LAPACKE_cgelsd( int matrix_order, lapack_int m, lapack_int n,
                           lapack_int nrhs, lapack_complex_float* a,
                           lapack_int lda, lapack_complex_float* b,
                           lapack_int ldb, float* s, float rcond,
                           lapack_int* rank )
{
    lapack_int info = 0;
    lapack_int lwork = -1;
    /* Additional scalars declarations for work arrays */
    lapack_int liwork;
    lapack_int lrwork;
    lapack_int* iwork = NULL;
    float* rwork = NULL;
    lapack_complex_float* work = NULL;
    lapack_int iwork_query;
    float rwork_query;
    lapack_complex_float work_query;
    if( matrix_order != LAPACK_COL_MAJOR && matrix_order != LAPACK_ROW_MAJOR ) {
        LAPACKE_xerbla( "LAPACKE_cgelsd", -1 );
        return -1;
    }
#ifndef LAPACK_DISABLE_NAN_CHECK
    /* Optionally check input matrices for NaNs */
    if( LAPACKE_cge_nancheck( matrix_order, m, n, a, lda ) ) {
        return -5;
    }
    if( LAPACKE_cge_nancheck( matrix_order, MAX(m,n), nrhs, b, ldb ) ) {
        return -7;
    }
    if( LAPACKE_s_nancheck( 1, &rcond, 1 ) ) {
        return -10;
    }
#endif
    /* Query optimal working array(s) size */
    info = LAPACKE_cgelsd_work( matrix_order, m, n, nrhs, a, lda, b, ldb, s,
                                rcond, rank, &work_query, lwork, &rwork_query,
                                &iwork_query );
    if( info != 0 ) {
        goto exit_level_0;
    }
    liwork = (lapack_int)iwork_query;
    lrwork = (lapack_int)rwork_query;
    lwork = LAPACK_C2INT( work_query );
    /* Allocate memory for work arrays */
    iwork = (lapack_int*)LAPACKE_malloc( sizeof(lapack_int) * liwork );
    if( iwork == NULL ) {
        info = LAPACK_WORK_MEMORY_ERROR;
        goto exit_level_0;
    }
    rwork = (float*)LAPACKE_malloc( sizeof(float) * lrwork );
    if( rwork == NULL ) {
        info = LAPACK_WORK_MEMORY_ERROR;
        goto exit_level_1;
    }
    work = (lapack_complex_float*)
        LAPACKE_malloc( sizeof(lapack_complex_float) * lwork );
    if( work == NULL ) {
        info = LAPACK_WORK_MEMORY_ERROR;
        goto exit_level_2;
    }
    /* Call middle-level interface */
    info = LAPACKE_cgelsd_work( matrix_order, m, n, nrhs, a, lda, b, ldb, s,
                                rcond, rank, work, lwork, rwork, iwork );
    /* Release memory and exit */
    LAPACKE_free( work );
exit_level_2:
    LAPACKE_free( rwork );
exit_level_1:
    LAPACKE_free( iwork );
exit_level_0:
    if( info == LAPACK_WORK_MEMORY_ERROR ) {
        LAPACKE_xerbla( "LAPACKE_cgelsd", info );
    }
    return info;
}

static int jiveL_initSDL(lua_State *L) {
	const SDL_VideoInfo *video_info;
	JiveSurface *srf, *splash, *icon;
	Uint16 splash_w, splash_h;

	/* logging */
	log_ui_draw = LOG_CATEGORY_GET("jivelite.ui.draw");
	log_ui = LOG_CATEGORY_GET("jivelite.ui");

	/* linux fbcon does not need a mouse */
	SDL_putenv("SDL_NOMOUSE=1");

	/* allow the screensaver */
	SDL_putenv("SDL_VIDEO_ALLOW_SCREENSAVER=1");

	/* initialise SDL */
	if (SDL_Init(SDL_INIT_VIDEO) < 0) {
		LOG_ERROR(log_ui_draw, "SDL_Init(V|T|A): %s\n", SDL_GetError());
		SDL_Quit();
		exit(-1);
	}

	/* report video info */
	if ((video_info = SDL_GetVideoInfo())) {
		LOG_INFO(log_ui_draw, "%d,%d %d bits/pixel %d bytes/pixel [R<<%d G<<%d B<<%d]", video_info->current_w, video_info->current_h, video_info->vfmt->BitsPerPixel, video_info->vfmt->BytesPerPixel, video_info->vfmt->Rshift, video_info->vfmt->Gshift, video_info->vfmt->Bshift);
		LOG_INFO(log_ui_draw, "Hardware acceleration %s available", video_info->hw_available?"is":"is not");
		LOG_INFO(log_ui_draw, "Window manager %s available", video_info->wm_available?"is":"is not");
	}

	/* Register callback for additional events (used for multimedia keys)*/
	SDL_EventState(SDL_SYSWMEVENT,SDL_ENABLE);
	SDL_SetEventFilter(filter_events);

	platform_init(L);

	/* open window */
	SDL_WM_SetCaption("JiveLite Alpha", "JiveLite Alpha");
	SDL_ShowCursor(SDL_DISABLE);
	SDL_EnableKeyRepeat (100, 100);
	SDL_EnableUNICODE(1);

	/* load the icon */
	icon = jive_surface_load_image("jive/app.png");
	if (icon) {
		jive_surface_set_wm_icon(icon);
		jive_surface_free(icon);
	}

	screen_w = video_info->current_w;
	screen_h = video_info->current_h;

	screen_bpp = video_info->vfmt->BitsPerPixel;

	splash = jive_surface_load_image("jive/splash.png");
	if (splash) {
		jive_surface_get_size(splash, &splash_w, &splash_h);
		if (video_info->wm_available) {
			screen_w = splash_w;
			screen_h = splash_h;
		}
	}

	srf = jive_surface_set_video_mode(screen_w, screen_h, screen_bpp, video_info->wm_available ? false : true);
	if (!srf) {
		SDL_Quit();
		exit(-1);
	}

	if (splash) {
		jive_surface_blit(splash, srf, MAX(0, (screen_w - splash_w) / 2), MAX(0, (screen_h - splash_h) / 2));
		jive_surface_flip(srf);
	}

	lua_getfield(L, 1, "screen");
	if (lua_isnil(L, -1)) {
		LOG_ERROR(log_ui_draw, "no screen table");

		SDL_Quit();
		exit(-1);
	}

	/* store screen surface */
	JiveSurface **p = (JiveSurface **)lua_newuserdata(L, sizeof(JiveSurface *));
	*p = jive_surface_ref(srf);
	luaL_getmetatable(L, "JiveSurface");
	lua_setmetatable(L, -2);

	lua_setfield(L, -2, "surface");

	lua_getfield(L, -1, "bounds");
	lua_pushinteger(L, screen_w);
	lua_rawseti(L, -2, 3);
	lua_pushinteger(L, screen_h);
	lua_rawseti(L, -2, 4);
	lua_pop(L, 2);

	/* background image */
	jive_background = jive_tile_fill_color(0x000000FF);

	/* jive.ui.style = {} */
//.........这里部分代码省略.........

// host stub function
void ops_par_loop_initialise_chunk_kernel_xx(char const *name, ops_block block,
                                             int dim, int *range, ops_arg arg0,
                                             ops_arg arg1) {

  // Timing
  double t1, t2, c1, c2;
  ops_arg args[2] = {arg0, arg1};

#ifdef CHECKPOINTING
  if (!ops_checkpointing_before(args, 2, range, 0))
    return;
#endif

  if (OPS_diags > 1) {
    ops_timing_realloc(0, "initialise_chunk_kernel_xx");
    OPS_kernels[0].count++;
    ops_timers_core(&c1, &t1);
  }

  // compute localy allocated range for the sub-block

  int start[2];
  int end[2];
#ifdef OPS_MPI
  sub_block_list sb = OPS_sub_block_list[block->index];
#endif // OPS_MPI

  int arg_idx[2];
  int arg_idx_base[2];
#ifdef OPS_MPI
  if (compute_ranges(args, 2, block, range, start, end, arg_idx) < 0)
    return;
#else // OPS_MPI
  for (int n = 0; n < 2; n++) {
    start[n] = range[2 * n];
    end[n] = range[2 * n + 1];
    arg_idx[n] = start[n];
  }
#endif
  for (int n = 0; n < 2; n++) {
    arg_idx_base[n] = arg_idx[n];
  }

  int dat0 = args[0].dat->elem_size;

  // set up initial pointers
  int base0 = args[0].dat->base_offset +
              (OPS_soa ? args[0].dat->type_size : args[0].dat->elem_size) *
                  start[0] * args[0].stencil->stride[0];
  base0 = base0 +
          (OPS_soa ? args[0].dat->type_size : args[0].dat->elem_size) *
              args[0].dat->size[0] * start[1] * args[0].stencil->stride[1];
#ifdef OPS_GPU
  int *p_a0 = (int *)((char *)args[0].data_d + base0);
#else
  int *p_a0 = (int *)((char *)args[0].data + base0);
#endif

  int *p_a1 = NULL;

  int x_size = MAX(0, end[0] - start[0]);
  int y_size = MAX(0, end[1] - start[1]);

  // initialize global variable with the dimension of dats
  xdim0 = args[0].dat->size[0];
  if (xdim0 != xdim0_initialise_chunk_kernel_xx_h) {
    xdim0_initialise_chunk_kernel_xx = xdim0;
    xdim0_initialise_chunk_kernel_xx_h = xdim0;
  }

// Halo Exchanges

#ifdef OPS_GPU
  ops_H_D_exchanges_device(args, 2);
#else
  ops_H_D_exchanges_host(args, 2);
#endif
  ops_halo_exchanges(args, 2, range);

#ifdef OPS_GPU
  ops_H_D_exchanges_device(args, 2);
#else
  ops_H_D_exchanges_host(args, 2);
#endif
  if (OPS_diags > 1) {
    ops_timers_core(&c2, &t2);
    OPS_kernels[0].mpi_time += t2 - t1;
  }

  initialise_chunk_kernel_xx_c_wrapper(p_a0, p_a1, arg_idx[0], arg_idx[1],
                                       x_size, y_size);

  if (OPS_diags > 1) {
    ops_timers_core(&c1, &t1);
    OPS_kernels[0].time += t1 - t2;
  }
#ifdef OPS_GPU
  ops_set_dirtybit_device(args, 2);
#else
//.........这里部分代码省略.........

void  MSXtank_mix2(int i, double vIn, double cIn[], long dt)
/*
**   Purpose: 2-compartment tank model                      
**
**   Input:   i = tank index
**            vIn = volume of inflow to tank (ft3)
**            cIn[] = concen. of inflow to tank
**            dt = current WQ time step (sec)
*/
{
    int     k, m, n;
    long    tstep,                     // Actual time step taken
            tstar;                     // Time to fill or drain a zone
    double  qIn,                       // Inflow rate
            qOut,                      // Outflow rate
            qNet;                      // Net flow rate
    double  c, c1, c2;                 // Species concentrations
    Pseg    seg1,                      // Mixing zone segment
            seg2;                      // Ambient zone segment

// --- find inflows & outflows 

    n = MSX.Tank[i].node;
    qNet = MSX.D[n];
    qIn = vIn/(double)dt;
    qOut = qIn - qNet;

// --- get segments for each zone

    k = MSX.Nobjects[LINK] + i;
    seg1 = MSX.FirstSeg[k];
    seg2 = MSX.LastSeg[k];

// --- case of no net volume change

    if ( fabs(qNet) < TINY ) return;

// --- case of net filling (qNet > 0)

    else if (qNet > 0.0)
    {

    // --- case where ambient zone empty & mixing zone filling 

        if (seg2->v <= 0.0)
        {
        // --- time to fill mixing zone
            tstar = (long) ((MSX.Tank[i].vMix - (seg1->v))/qNet);
            tstep = MIN(dt, tstar);

            for (m=1; m<=MSX.Nobjects[SPECIES]; m++)
            {
                if ( MSX.Species[m].type != BULK ) continue;

            // --- new quality in mixing zone
                c = seg1->c[m]; 
                if (seg1->v > 0.0) seg1->c[m] += qIn*tstep*(cIn[m]-c)/(seg1->v);
                else seg1->c[m] = cIn[m];
                seg1->c[m] = MAX(0.0, seg1->c[m]);
                seg2->c[m] = 0.0;
            }

        // --- new volume of mixing zone
            seg1->v += qNet*tstep;

        // --- time during which ambient zone fills
            dt -= tstep;
        }

    // --- case where mixing zone full & ambient zone filling 
        if (dt > 1)
        {
            for (m=1; m<=MSX.Nobjects[SPECIES]; m++)
            {
                if ( MSX.Species[m].type != BULK ) continue;

            // --- new quality in mixing zone
                c1 = seg1->c[m]; 
                seg1->c[m] += qIn * dt * (cIn[m] - c1) / (seg1->v);
                seg1->c[m] = MAX(0.0, seg1->c[m]);

            // --- new quality in ambient zone
		        c2 = seg2->c[m];
                if (seg2->v <= 0.0)
                    seg2->c[m] = seg1->c[m];
                else
                    seg2->c[m] += qNet * dt * ((seg1->c[m]) - c2) / (seg2->v);
                seg2->c[m] = MAX(0.0, seg2->c[m]);
	        }

        // --- new volume of ambient zone
            seg2->v += qNet*dt;
        }
        if ( seg1->v > 0.0 ) MSXchem_equil(NODE, seg1->c);
        if ( seg2->v > 0.0 ) MSXchem_equil(NODE, seg2->c);
    }

// --- case of net emptying (qnet < 0)

    else if ( qNet < 0.0 && seg1->v > 0.0 )
//.........这里部分代码省略.........

static int mag_manacost(struct char_data *ch, int spellnum)
{
  return MAX(SINFO.mana_max - (SINFO.mana_change *
		    (GET_LEVEL(ch) - SINFO.min_level[(int) GET_CLASS(ch)])),
	     SINFO.mana_min);
}

static void
mx_texture_frame_paint_texture_internal (CoglHandle  material,
                                         CoglHandle  texture,
                                         guint8      opacity,
                                         gfloat      top,
                                         gfloat      right,
                                         gfloat      bottom,
                                         gfloat      left,
                                         gfloat      width,
                                         gfloat      height)
{
  gfloat tex_width, tex_height;
  gfloat ex, ey;
  gfloat tx1, ty1, tx2, ty2;

  /* apply opacity */
  cogl_material_set_color4ub (material, opacity, opacity, opacity, opacity);

  /* add the texture */
  cogl_material_set_layer (material, 0, texture);

  /* set the source */
  cogl_set_source (material);

  tex_width  = cogl_texture_get_width (texture);
  tex_height = cogl_texture_get_height (texture);

  /* simple stretch */
  if (left == 0 && right == 0 && top == 0
      && bottom == 0)
    {
      cogl_rectangle (0, 0, width, height);
      return;
    }

  tx1 = left / tex_width;
  tx2 = (tex_width - right) / tex_width;
  ty1 = top / tex_height;
  ty2 = (tex_height - bottom) / tex_height;

  ex = width - right;
  if (ex < left)
    ex = left;

  ey = height - bottom;
  if (ey < top)
    ey = top;


  {
    float rectangles[] =
    {
      /* top left corner */
      0, 0,
      left, top,
      0.0, 0.0,
      tx1, ty1,

      /* top middle */
      left, 0,
      MAX (left, ex), top,
      tx1, 0.0,
      tx2, ty1,

      /* top right */
      ex, 0,
      MAX (ex + right, width), top,
      tx2, 0.0,
      1.0, ty1,

      /* mid left */
      0, top,
      left,  ey,
      0.0, ty1,
      tx1, ty2,

      /* center */
      left, top,
      ex, ey,
      tx1, ty1,
      tx2, ty2,

      /* mid right */
      ex, top,
      MAX (ex + right, width), ey,
      tx2, ty1,
      1.0, ty2,

      /* bottom left */
      0, ey,
      left, MAX (ey + bottom, height),
      0.0, ty2,
      tx1, 1.0,

      /* bottom center */
      left, ey,
      ex, MAX (ey + bottom, height),
      tx1, ty2,
      tx2, 1.0,

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

void Blob2D::resize(short destLeft, short destTop, short destWidth, short destHeight)
{
	if( width == 0 && height == 0)
	{
		size_t siz = Bitmap::size(destWidth,destHeight);
		if( alloc_size < siz )
		{
			if(ptr)
				mem_del(ptr);
			ptr = (Bitmap *)mem_add((alloc_size = siz + DEFAULT_BLOB2D_INC));
		}
		memset(ptr, BACKGROUND_COLOR, siz);
		ptr->init(destWidth, destHeight);
	}
	else if( destWidth == 0 && destHeight == 0 )
	{
		err_when(!ptr);
		ptr->init(destWidth, destHeight);
	}
	else if( left_edge == destLeft && top_edge == destTop && width == destWidth )
	{
		if( destHeight == height )
			return;		// unchange

		size_t siz = Bitmap::size(destWidth,destHeight);
		if( alloc_size < siz )
		{
			ptr = (Bitmap *)mem_resize(ptr, (alloc_size = siz + DEFAULT_BLOB2D_INC));
		}

		if( destHeight > height && destWidth > 0 )
		{
			int y2 = top_edge+height;	// must keep the old instant
			height = destHeight;			// as height must be change to make fill_area correct

			// extend, then fill the new with the background color
			fill_area(destLeft, top_edge+height, destLeft+destWidth-1, destTop+destHeight-1, BACKGROUND_COLOR, 0);
		}
	}
	else if( left_edge <= destLeft && top_edge <= destTop &&
		left_edge+width >= destLeft+destWidth && top_edge+height >= destTop+destHeight)
	{
		// clipping
		unsigned char *src = ptr->get_ptr(destLeft-left_edge, destTop-top_edge);
		int srcPitch = ptr->get_pitch();
		ptr->init(destWidth, destHeight);
		unsigned char *dest = ptr->get_ptr();
		int destPitch = ptr->get_pitch();

		for(int y = 0; y < destHeight; ++y, src += srcPitch, dest += destPitch )
			memmove(dest, src, destWidth);

		// ptr = (Bitmap *)mem_resize(ptr, ptr->size());
	}
	else
	{
		// general resize, new another buffer
		// copy range, intersection of two area :
		short copyLeft, copyTop, copyWidth, copyHeight;
		copyLeft = MAX(destLeft, left_edge);
		copyTop = MAX(destTop, top_edge);
		copyWidth = MIN(destLeft + destWidth, left_edge + width) - copyLeft;
		copyHeight = MIN(destTop + destHeight, top_edge + height) - copyTop;

		{
			size_t siz = Bitmap::size(destWidth, destHeight);
			Bitmap *newPtr = (Bitmap *)mem_add(siz + DEFAULT_BLOB2D_INC);
			memset(newPtr, BACKGROUND_COLOR, siz);
			newPtr->init(destWidth, destHeight);

			if( copyWidth > 0 && copyHeight > 0 )
			{
				int yCount = 0;
				unsigned char *src = ptr->get_ptr(copyLeft-left_edge, yCount+copyTop-top_edge );
				unsigned char *dest = newPtr->get_ptr(copyLeft-destLeft, yCount+copyTop-destTop );

				for( ; yCount < copyHeight; ++yCount, src += ptr->get_pitch(), dest += ptr->get_pitch() )
				{
					// unsigned char *src = (yCount+copyTop-top_edge)*width + copyLeft-left_edge;
					// unsigned char *dest = (yCount+copyTop-destTop)*destWdith + copyLeft-destLeft;
					memcpy(dest, src, copyWidth);
				}
			}

			// assign to the newPtr now
			left_edge = destLeft;
			top_edge = destTop;
			width = destWidth;
			height = destHeight;
			if(ptr)
				mem_del(ptr);
			ptr = newPtr;
		}

		// fill rest area with background color
		if( top_edge < copyTop && width > 0)
		{
			fill_area(left_edge, top_edge, left_edge+width-1, copyTop, BACKGROUND_COLOR, 0 );
		}

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

STATIC
UINT64*
GetBlockEntryListFromAddress (
  IN  UINT64       *RootTable,
  IN  UINT64        RegionStart,
  OUT UINTN        *TableLevel,
  IN OUT UINT64    *BlockEntrySize,
  OUT UINT64      **LastBlockEntry
  )
{
  UINTN   RootTableLevel;
  UINTN   RootTableEntryCount;
  UINT64 *TranslationTable;
  UINT64 *BlockEntry;
  UINT64 *SubTableBlockEntry;
  UINT64  BlockEntryAddress;
  UINTN   BaseAddressAlignment;
  UINTN   PageLevel;
  UINTN   Index;
  UINTN   IndexLevel;
  UINTN   T0SZ;
  UINT64  Attributes;
  UINT64  TableAttributes;

  // Initialize variable
  BlockEntry = NULL;

  // Ensure the parameters are valid
  if (!(TableLevel && BlockEntrySize && LastBlockEntry)) {
    ASSERT_EFI_ERROR (EFI_INVALID_PARAMETER);
    return NULL;
  }

  // Ensure the Region is aligned on 4KB boundary
  if ((RegionStart & (SIZE_4KB - 1)) != 0) {
    ASSERT_EFI_ERROR (EFI_INVALID_PARAMETER);
    return NULL;
  }

  // Ensure the required size is aligned on 4KB boundary and not 0
  if ((*BlockEntrySize & (SIZE_4KB - 1)) != 0 || *BlockEntrySize == 0) {
    ASSERT_EFI_ERROR (EFI_INVALID_PARAMETER);
    return NULL;
  }

  T0SZ = ArmGetTCR () & TCR_T0SZ_MASK;
  // Get the Table info from T0SZ
  GetRootTranslationTableInfo (T0SZ, &RootTableLevel, &RootTableEntryCount);

  // If the start address is 0x0 then we use the size of the region to identify the alignment
  if (RegionStart == 0) {
    // Identify the highest possible alignment for the Region Size
    BaseAddressAlignment = LowBitSet64 (*BlockEntrySize);
  } else {
    // Identify the highest possible alignment for the Base Address
    BaseAddressAlignment = LowBitSet64 (RegionStart);
  }

  // Identify the Page Level the RegionStart must belong to. Note that PageLevel
  // should be at least 1 since block translations are not supported at level 0
  PageLevel = MAX (3 - ((BaseAddressAlignment - 12) / 9), 1);

  // If the required size is smaller than the current block size then we need to go to the page below.
  // The PageLevel was calculated on the Base Address alignment but did not take in account the alignment
  // of the allocation size
  while (*BlockEntrySize < TT_BLOCK_ENTRY_SIZE_AT_LEVEL (PageLevel)) {
    // It does not fit so we need to go a page level above
    PageLevel++;
  }

  //
  // Get the Table Descriptor for the corresponding PageLevel. We need to decompose RegionStart to get appropriate entries
  //

  TranslationTable = RootTable;
  for (IndexLevel = RootTableLevel; IndexLevel <= PageLevel; IndexLevel++) {
    BlockEntry = (UINT64*)TT_GET_ENTRY_FOR_ADDRESS (TranslationTable, IndexLevel, RegionStart);

    if ((IndexLevel != 3) && ((*BlockEntry & TT_TYPE_MASK) == TT_TYPE_TABLE_ENTRY)) {
      // Go to the next table
      TranslationTable = (UINT64*)(*BlockEntry & TT_ADDRESS_MASK_DESCRIPTION_TABLE);

      // If we are at the last level then update the last level to next level
      if (IndexLevel == PageLevel) {
        // Enter the next level
        PageLevel++;
      }
    } else if ((*BlockEntry & TT_TYPE_MASK) == TT_TYPE_BLOCK_ENTRY) {
      // If we are not at the last level then we need to split this BlockEntry
      if (IndexLevel != PageLevel) {
        // Retrieve the attributes from the block entry
        Attributes = *BlockEntry & TT_ATTRIBUTES_MASK;

        // Convert the block entry attributes into Table descriptor attributes
        TableAttributes = TT_TABLE_AP_NO_PERMISSION;
        if (Attributes & TT_NS) {
          TableAttributes = TT_TABLE_NS;
        }

        // Get the address corresponding at this entry
//.........这里部分代码省略.........

static void loop()
{
    static const uint8_t compass_count = compass.get_count();
    static float min[COMPASS_MAX_INSTANCES][3];
    static float max[COMPASS_MAX_INSTANCES][3];
    static float offset[COMPASS_MAX_INSTANCES][3];

    if ((AP_HAL::micros() - timer) > 100000L) {
        timer = AP_HAL::micros();
        compass.read();
        const uint32_t read_time = AP_HAL::micros() - timer;

        for (uint8_t i = 0; i < compass_count; i++) {
            float heading;

            hal.console->printf("Compass #%u: ", i);

            if (!compass.healthy()) {
                hal.console->printf("not healthy\n");
                continue;
            }

            Matrix3f dcm_matrix;
            // use roll = 0, pitch = 0 for this example
            dcm_matrix.from_euler(0, 0, 0);
            heading = compass.calculate_heading(dcm_matrix, i);

            const Vector3f &mag = compass.get_field(i);

            // capture min
            min[i][0] = MIN(mag.x, min[i][0]);
            min[i][1] = MIN(mag.y, min[i][1]);
            min[i][2] = MIN(mag.z, min[i][2]);

            // capture max
            max[i][0] = MAX(mag.x, max[i][0]);
            max[i][1] = MAX(mag.y, max[i][1]);
            max[i][2] = MAX(mag.z, max[i][2]);

            // calculate offsets
            offset[i][0] = -(max[i][0] + min[i][0]) / 2;
            offset[i][1] = -(max[i][1] + min[i][1]) / 2;
            offset[i][2] = -(max[i][2] + min[i][2]) / 2;

            // display all to user
            hal.console->printf("Heading: %.2f (%3d, %3d, %3d)",
                                (double)ToDeg(heading),
                                (int)mag.x,
                                (int)mag.y,
                                (int)mag.z);

            // display offsets
            hal.console->printf(" offsets(%.2f, %.2f, %.2f)",
                                (double)offset[i][0],
                                (double)offset[i][1],
                                (double)offset[i][2]);

            hal.console->printf(" t=%u", (unsigned)read_time);

            hal.console->printf("\n");
        }
    } else {
        hal.scheduler->delay(1);
    }
}

static void fmtint(char *buffer, size_t *currlen, size_t maxlen,
		    long value, int base, int min, int max, int flags)
{
	int signvalue = 0;
	unsigned long uvalue;
	char convert[20];
	int place = 0;
	int spadlen = 0; /* amount to space pad */
	int zpadlen = 0; /* amount to zero pad */
	int caps = 0;
	
	if (max < 0)
		max = 0;
	
	uvalue = value;
	
	if(!(flags & DP_F_UNSIGNED)) {
		if( value < 0 ) {
			signvalue = '-';
			uvalue = -value;
		} else {
			if (flags & DP_F_PLUS)  /* Do a sign (+/i) */
				signvalue = '+';
			else if (flags & DP_F_SPACE)
				signvalue = ' ';
		}
	}
  
	if (flags & DP_F_UP) caps = 1; /* Should characters be upper case? */

	do {
		convert[place++] =
			(caps? "0123456789ABCDEF":"0123456789abcdef")
			[uvalue % (unsigned)base  ];
		uvalue = (uvalue / (unsigned)base );
	} while(uvalue && (place < 20));
	if (place == 20) place--;
	convert[place] = 0;

	zpadlen = max - place;
	spadlen = min - MAX (max, place) - (signvalue ? 1 : 0);
	if (zpadlen < 0) zpadlen = 0;
	if (spadlen < 0) spadlen = 0;
	if (flags & DP_F_ZERO) {
		zpadlen = MAX(zpadlen, spadlen);
		spadlen = 0;
	}
	if (flags & DP_F_MINUS) 
		spadlen = -spadlen; /* Left Justifty */

#ifdef DEBUG_SNPRINTF
	printf("zpad: %d, spad: %d, min: %d, max: %d, place: %d\n",
	       zpadlen, spadlen, min, max, place);
#endif

	/* Spaces */
	while (spadlen > 0) {
		dopr_outch (buffer, currlen, maxlen, ' ');
		--spadlen;
	}

	/* Sign */
	if (signvalue) 
		dopr_outch (buffer, currlen, maxlen, signvalue);

	/* Zeros */
	if (zpadlen > 0) {
		while (zpadlen > 0) {
			dopr_outch (buffer, currlen, maxlen, '0');
			--zpadlen;
		}
	}

	/* Digits */
	while (place > 0) 
		dopr_outch (buffer, currlen, maxlen, convert[--place]);
  
	/* Left Justified spaces */
	while (spadlen < 0) {
		dopr_outch (buffer, currlen, maxlen, ' ');
		++spadlen;
	}
}

/******************************************************************************
Description.: this is the main worker thread
              it loops forever, grabs a fresh frame and calculates focus
Input Value.:
Return Value:
******************************************************************************/
void *worker_thread( void *arg ) {
  int frame_size=0;
  double sv=-1.0, max_sv=100.0, delta=500;
  int focus=255, step=10, max_focus=100, search_focus=1;

  if ( (frame = malloc(256*1024)) == NULL ) {
    OPRINT("not enough memory for worker thread\n");
    exit(EXIT_FAILURE);
  }

  /* set cleanup handler to cleanup allocated ressources */
  pthread_cleanup_push(worker_cleanup, NULL);

  while ( !pglobal->stop ) {
    DBG("waiting for fresh frame\n");
    pthread_cond_wait(&pglobal->db_update, &pglobal->db);

    /* read buffer */
    frame_size = pglobal->size;
    memcpy(frame, pglobal->buf, frame_size);

    pthread_mutex_unlock( &pglobal->db );

    /* process frame */
    sv = getFrameSharpnessValue(frame, frame_size);
    DBG("sharpness is: %f\n", sv);

    if ( search_focus || (ABS(sv-max_sv) > delta) ) {
      DBG("adjusting focus: %d\n", focus);

      /* entered because focus changed */
      if ( !search_focus ) {
        DBG("starting to search for focus\n");
        max_focus    = 255;
        focus        = 255;
        max_sv       = -1.0;
        search_focus = 1;
      }

      if ( focus <= 0 ) {
        focus = max_focus;
        DBG("max focus found at: %d\n", max_focus);
        search_focus = 0;
      }

      if ( search_focus ) {
        if (sv > max_sv) {
          /* sharpness is better then max now */
          DBG("found better focus at: %d\n", focus);
          max_focus = focus;
          max_sv = sv;
        }

        focus = MIN(MAX(focus-step,0), 255);
        DBG("decrement focus now to: %d\n", focus);
        //focus = pglobal->in.cmd(IN_CMD_FOCUS_SET, focus);
      }
    }

    if ( (delay > 0) && !search_focus ) {
      usleep(1000*delay);
    }
  }

  pthread_cleanup_pop(1);

  return NULL;
}

//.........这里部分代码省略.........
                    /* general iteration */

                    vbs = model->MOS3type * ( 
                        *(ckt->CKTrhsOld+here->MOS3bNode) -
                        *(ckt->CKTrhsOld+here->MOS3sNodePrime));
                    vgs = model->MOS3type * ( 
                        *(ckt->CKTrhsOld+here->MOS3gNode) -
                        *(ckt->CKTrhsOld+here->MOS3sNodePrime));
                    vds = model->MOS3type * ( 
                        *(ckt->CKTrhsOld+here->MOS3dNodePrime) -
                        *(ckt->CKTrhsOld+here->MOS3sNodePrime));

                /* now some common crunching for some more useful quantities */



            /*
             * now all the preliminaries are over - we can start doing the
             * real work
             */
            vbd = vbs - vds;
            vgd = vgs - vds;
            vgb = vgs - vbs;

            /* bulk-source and bulk-drain doides
             * here we just evaluate the ideal diode current and the
             * correspoinding derivative (conductance).
             */
	    if(vbs <= 0) {
                lgbs = SourceSatCur/vt;
                lgbs += ckt->CKTgmin;
		lgbs2 = lgbs3 = 0;
            } else {
                evbs = exp(MIN(MAX_EXP_ARG,vbs/vt));
                lgbs = SourceSatCur*evbs/vt + ckt->CKTgmin;
		lgbs2 = model->MOS3type *0.5 * (lgbs - ckt->CKTgmin)/vt;
		lgbs3 = model->MOS3type *lgbs2/(vt*3);

            }
            if(vbd <= 0) {
                lgbd = DrainSatCur/vt;
                lgbd += ckt->CKTgmin;
		lgbd2 = lgbd3 = 0;
            } else {
                ebd = exp(MIN(MAX_EXP_ARG,vbd/vt));
                lgbd = DrainSatCur*ebd/vt +ckt->CKTgmin;
		lgbd2 = model->MOS3type *0.5 * (lgbd - ckt->CKTgmin)/vt;
		lgbd3 = model->MOS3type *lgbd2/(vt*3);
            }


            /* now to determine whether the user was able to correctly
             * identify the source and drain of his device
             */
            if(vds >= 0) {
                /* normal mode */
                here->MOS3mode = 1;
            } else {
                /* inverse mode */
                here->MOS3mode = -1;
            }

            {
            /*
             * subroutine moseq3(vds,vbs,vgs,gm,gds,gmbs,
             * qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)

//.........这里部分代码省略.........
     Coditional Reduction Functions
  */
  D