stats_utils.c

卫星仿真软件 卫星仿真软件 卫星仿真软件

	}
	for(j = left[0]; j < right_edge[0]; j++) {
	  increment_intensity(row,j,g);
	}
      }
      prev_row = row;
    }
  }

  /*
   * Footprint caps one of the poles if projection of center vector c onto
   * x,y plane has magnitude less than magnitude of projection of circle basis
   * vector u (u always points in direction of north pole from c).
   */
  if ( pc->x*pc->x+pc->y*pc->y < pu->x*pu->x+pu->y*pu->y ) {
    cap_pole( pc, pu, pv, current_projection, g);
  }
}

/*
 * Increment the intensity image at given row and column.
 * (row ranges from 0 to g->height-1,
 *  column from 0 to g->width-1.)
 */
static inline void
increment_intensity(int row, int column, grid *g)
{
  unsigned int k;

  k = g->width * row + column;

  if ( 0 == (g->data[k])++ ) {
    (g->noaccess[k]) = 0;
    (g->covered[row])++;
    (g->count)++;
  }
}

/*
 * For time t, return the spherical coordinates of the point
 * c+cos(t)u+sin(t)v on circle with center c and basis vectors u and v.
 */
static void
circle_point(SphericalCoordinates *point, CartesianCoordinates *pc,
	     CartesianCoordinates *pu, CartesianCoordinates *pv, double t )
{
  CartesianCoordinates temp;
  float ct = cos(t), st = sin(t);

  temp.x = pc->x+ct*pu->x+st*pv->x;
  temp.y = pc->y+ct*pu->y+st*pv->y;
  temp.z = pc->z+ct*pu->z+st*pv->z;
  cartesian_to_spherical(point, &temp);
}

/*
 * Given point on unit sphere in spherical coordinates, return corresponding
 * index into coverage intensity grid.
 *        grid_index[0]=column index, grid_index[1]=row index.
 */
static void
intensity_index(int grid_index[2], SphericalCoordinates *point,
		int current_projection, grid *g)
{
  double proj[2];

  switch (current_projection) {
  case UNPROJECTED:
    project_unprojected(proj, point);
    break;
  case SINUSOIDAL:
    project_sinusoidal(proj, point);
    break;
  case SPHERICAL:
    project_spherical(proj,point);
    break;
  case CYLINDRICAL:
  default:
    project_cylindrical(proj, point);
  }

  /* It's assumed that proj[0] ranges from -PI to PI no matter what. */
  grid_index[0]=((int) ((proj[0]+PI)*(g->width)/TWOPI));
  grid_index[1]=((int) (proj[1]*(g->height-1)/PI));
}

/*
 * Increment intensity image for a polar cap.
 *
 * The routine increments those image pixels corresponding to the largest
 * disk centered at the pole contained within the footprint.  We're helped
 * here by the fact that u, when based at the tip of the center vector c,
 * always points in the direction of the north pole.  So for north pole
 * caps, we simply increment rows 0 through the row corresponding to c+u.
 *
 * For south poles, the vector c-u is the closest point in the footprint
 * to the pole, so we increment that row through the final row.
 */
static void
cap_pole(CartesianCoordinates *pc, CartesianCoordinates *pu,
	 CartesianCoordinates *pv, int current_projection, grid *g )
{
  SphericalCoordinates point;
  int grid_index[2], left[2], right[2], bottom, i, j;
  const double dphi=PI/g->height;

  if (pc->z > 0.0) {  /* Cap covers north pole. */
    circle_point(&point, pc, pu, pv, 0.0);
    intensity_index(grid_index, &point, current_projection, g);

    if (debug) {
      fprintf(stderr, "North pole cap: begin %d, end %d\n", 0, grid_index[1]);
      fprintf(stderr, "Size of cap: theta=%f.2\n", point.theta);
    }

    point.phi = 0.0;

    for (i = 0; i < grid_index[1]+1; i++) {

      if (debug) fprintf(stderr, "row %d  ", i);

      intensity_edges(left, right, &point, current_projection, g);
      if (i == 0 && right[0] <= left[0]) right[0] = left[0] + 1;

      if (debug) fprintf(stderr, "columns %d to %d\n", left[0], right[0]);

      for (j = left[0]; j < right[0]; j++) {
	increment_intensity(i,j,g);
      }
      point.phi += dphi;
    }
  } else { /* South pole. But u still points towards north pole. */
    circle_point(&point, pc, pu, pv, PI); /* t=PI gets point on other side */
    intensity_index(grid_index, &point, current_projection, g); /* of south pole. */
    bottom = g->height;

    if (debug) fprintf(stderr, "South pole cap: begin %d, end %d\n",
		       grid_index[1]+1, bottom);

    for (i = grid_index[1]+1; i< bottom; i++) {

      if (debug) fprintf(stderr, "row %d  ", i);

      intensity_edges(left, right, &point, current_projection, g);
      if (i == bottom - 1 && right[0] <= left[0]) right[0] = left[0] + 1;

      if (debug) fprintf(stderr, "columns %d to %d\n", left[0], right[0]);

      for(j = left[0]; j < right[0]; j++) {
	increment_intensity(i,j,g);
      }
      point.phi += dphi;
    }
  }
}

/*
 * Given spherical coordinates on unit sphere, compute area-preserving
 * sinusoidal 2-d projection. Map is (1,theta,phi) -> (theta*sin(phi), phi).
 * proj[0] in [-Pi, Pi], proj[1] in [0, Pi].
 */
static void
project_sinusoidal(double proj[2], SphericalCoordinates *point)
{
  proj[0]=point->theta*sin(point->phi);
  proj[1]=point->phi;
}

/*
 * Given spherical coordinates on unit sphere, compute "area-preserving"
 * rectangular 2-d projection. Map is (1,theta,phi) -> (theta, Pi/2 * (1-cos(phi))).
 * proj[0] in [-Pi, Pi], proj[1] in [0, Pi] with phi=0 corresponding to
 * y=0.  (dA on sphere = 2/Pi dx dy)
 */
static void
project_cylindrical(double proj[2], SphericalCoordinates *point)
{
  proj[0]=point->theta;
  proj[1]=HALFPI*(1.0-cos(point->phi));
}

/*
 * Given spherical coordinates on unit sphere, compute unprojected
 * rectangular 2-d projection. Map is (1,theta,phi) -> (theta, phi).
 * proj[0] in [-Pi, Pi], proj[1] in [0, Pi].
 */
static void
project_unprojected(double proj[2], SphericalCoordinates *point)
{
  proj[0]=point->theta;
  proj[1]=point->phi;
}

/*
 * Given spherical coordinates on unit sphere, compute spherical
 * globes, two to a rectangle. Map is (1,theta,phi) -> (theta, phi).
 * proj[0] in [-Pi, Pi], proj[1] in [0, Pi].
 */
static void
project_spherical(double proj[2], SphericalCoordinates *point)
{
  if (point->theta < 0) {
    proj[0] = -HALFPI + HALFPI * sin(point->theta + HALFPI) * cos(HALFPI-point->phi);
  } else {
    proj[0] = HALFPI - HALFPI * sin(point->theta + HALFPI) * cos(HALFPI-point->phi);
  }
  proj[1]=HALFPI*(1.0-cos(point->phi));
}


/*
 * Given lat/lon coordinates on sphere, compute area-preserving
 * sinusoidal 2-d projection. Map is (lon, lat) -> (lon*sin(90-lat), 90-lat).
 * proj[0] in [-180, 180], proj[1] in [0,180].
 */
void
project_latlon_sinusoidal(double proj[2], LatLon *point)
{
  proj[0]=point->lon*(cos(DEG_TO_RAD*(point->lat)));
  proj[1]=90.0-point->lat;
}

/*
 * Given lat/lon coordinates on sphere, compute "area-preserving"
 * rectangular 2-d projection. Map is
 * (lon, lat) -> (lon, 90(1-cos(90-lat))).
 * proj[0] in [-180,180], proj[1] in [0,180] with lon=0 corresponding to
 * y=0.  (dA on sphere = 2/Pi dx dy)
 */
void
project_latlon_cylindrical(double proj[2], LatLon *point)
{
  proj[0]=point->lon;
  proj[1]=90.0*(1-sin(DEG_TO_RAD*(point->lat)));
}

/*
 * Given lat/lon coordinates on sphere, compute rectangular
 * unprojected 2-d projection. Map is (lon, lat) -> (lon*sin(90-lat), 90-lat).
 * proj[0] in [-180, 180], proj[1] in [0,180].
 */
void
project_latlon_unprojected(double proj[2], LatLon *point)
{
  proj[0]=point->lon;
  proj[1]=90.0-point->lat;
}

/*
 * Given lat/lon coordinates on sphere, compute spherical globes
 * on 2-d projection. Map is (lon, lat) -> broken into two hemispheres.
 * proj[0] in [-180, 180], proj[1] in [0,180].
 */
void
project_latlon_spherical(double proj[2], LatLon *point)
{
  if (point->lon < 0) {
    proj[0] =-90 + 90.0*sin(DEG_TO_RAD*(point->lon + 90)) * cos(DEG_TO_RAD*(point->lat));
  } else {
    proj[0] = 90 - 90.0*sin(DEG_TO_RAD*(point->lon + 90)) * cos(DEG_TO_RAD*(point->lat));
  }
  proj[1] = 90.0*(1-sin(DEG_TO_RAD*(point->lat)));
}

/*
 * Compute the left and right edge points in the projection map at the same
 * height as the projection of a point given in spherical coords (r, theta,
 * phi).
 *
 * left, right[0] = x (row) position
 * left, right[1] = y (column) position
 */
void
intensity_edges(int left[2], int right[2], SphericalCoordinates *point,
		int current_projection, grid *g)
{
  SphericalCoordinates temp;

  temp.r = point->r;
  temp.theta = PI;
  temp.phi = point->phi;
  intensity_index(right, &temp, current_projection, g);
  left[1]=right[1];
  left[0]=g->width-right[0];
}