stats_utils.c

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

/*
 *****************************************************
 *
 *  SaVi by Robert Thurman ([e-mail]) and
 *          Patrick Worfolk ([e-mail]).
 *
 *  Copyright (c) 1997 by The Geometry Center.
 *  This file is part of SaVi.  SaVi is free software;
 *  you can redistribute it and/or modify it only under
 *  the terms given in the file COPYRIGHT which you should
 *  have received along with this file.  SaVi may be
 *  obtained from:
 *  http://savi.sourceforge.net/
 *  http://www.geom.uiuc.edu/locate/SaVi
 *
 *****************************************************
 *
 * stats_utils.c
 *
 * $Id: stats_utils.c,v 1.18 2005/02/12 15:52:33 lloydwood Exp $
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "sats.h"
#include "constants.h"
#include "globals.h"
#include "stats_utils.h"
#include "orbit_utils.h"
#include "utils.h"


static void orient_circle(double t, CartesianCoordinates *pc,
			  CartesianCoordinates *pu,
			  CartesianCoordinates *pv,
			  const CentralBody *pcb);
static void footprint_circle(CartesianCoordinates *pc,
			     CartesianCoordinates *pu,
			     CartesianCoordinates *pv,
			     int fp_angle_type, double footprint_angle,
			     CartesianCoordinates *py, const CentralBody *pcb);
static double footprint_circle_radius(double l,
				      int fp_angle_type,
				      double footprint_angle);
static void intensity_circle_footprint(CartesianCoordinates *pc,
				       CartesianCoordinates *pu,
				       CartesianCoordinates *pv,
				       int current_projection, grid *g);

static inline void increment_intensity (int row, int column, grid *g);
static void circle_point(SphericalCoordinates *point,
			 CartesianCoordinates *pc,
			 CartesianCoordinates *pu,
			 CartesianCoordinates *pv,
			 double t);
static void intensity_index(int grid_index[2], SphericalCoordinates *point,
			    int current_projection, grid *g);
static void cap_pole(CartesianCoordinates *pc, CartesianCoordinates *pu,
		     CartesianCoordinates *pv,
		     int current_projection, grid *g);
static void project_sinusoidal(double proj[2], SphericalCoordinates *point);
static void project_cylindrical(double proj[2], SphericalCoordinates *point);
static void project_unprojected(double proj[2], SphericalCoordinates *point);
static void project_spherical(double proj[2], SphericalCoordinates *point);

/*
 * Initialize the coverage grid
 */
grid*
create_grid(int h, int w)
{
  grid *pg;

  pg=(grid *) malloc(sizeof(grid));

  pg->data=(int *) malloc(h * w *sizeof(int));
  pg->noaccess=(int *) malloc(h * w * sizeof(int));
  pg->covered=(int *) malloc(h * sizeof(int));

  pg->height=h;
  pg->width=w;
  pg->count = 0;

  return pg;
}

/*
 * Free space allocated to a coverage grid.
 */
void
destroy_grid(grid *g)
{
  free(g->data);
  free(g->noaccess);
  free(g->covered);

  free(g);
}

/*
 * Fill out the coverage grid.
 */
void
fill_grid(const Satellite_list SL, int current_projection, int fp_angle_type,
	   double footprint_angle, const CentralBody *pcb, grid *g)
{
  CartesianCoordinates c, u, v;
  Satellite_list sl;

  for (sl = SL; NULL != sl; sl = sl->next) {
    footprint_circle(&c, &u, &v, fp_angle_type, footprint_angle,
		     &(sl->s->x_C), pcb);
    orient_circle(sl->s->t, &c, &u, &v, pcb);
    intensity_circle_footprint(&c, &u, &v, current_projection, g);
  }
}

/*
 * Given orbital elements, compute circular footprint of a satellite.
 *
 * The circle is normalized to lie on a unit radius reference sphere, and
 * is specified by a center c and basis vectors u and v, in normalized
 * geocentric coordinates.  The parametrization is thus c+cos(t)u + sin(t)v.
 */
static void
footprint_circle(CartesianCoordinates *pc, CartesianCoordinates *pu,
		 CartesianCoordinates *pv,
		 int fp_angle_type, double footprint_angle,
		 CartesianCoordinates *py, const CentralBody *pcb)
{
  SphericalCoordinates z;
  CartesianCoordinates w;
  double ly, l, d, r, ip;

/* The center first */

  ly = norm(py);
  l = ly/pcb->radius;
  r = footprint_circle_radius( l, fp_angle_type, footprint_angle );
  d = sqrt(1-r*r);
  ip = d / ly;
  pc->x = py->x*ip;
  pc->y = py->y*ip;
  pc->z = py->z*ip;

/* Now the "up" basis vector. */

  cartesian_to_spherical(&z, py);
  z.r = 1;
  z.phi -= asin(r);

  if (debug) fprintf(stderr, "r=%.4f dphi=%.4f y=(%.2f, %.2f, %.2f)\n", \
	  r, asin(r), z.r, z.theta, z.phi);

  spherical_to_cartesian(&w, &z);
  pu->x = w.x-pc->x;
  pu->y = w.y-pc->y;
  pu->z = w.z-pc->z;

/* The final basis vector we can take to be the cross-product c x u,
 * normalized to have same length as u (that is, r).  Notice, then, that
 * v points to the "left". */

  cross_product( pv, pc, pu );
  pv->x /= d;
  pv->y /= d;
  pv->z /= d;
}

/*
 * orient_circle
 *
 * For displaying onto the coverage map, the footprint circle vectors need
 * to be rotated around the z-axis so that Longitude_Center_Line appears
 * at geocentric equatorial theta=0.
 */
static void
orient_circle(double t, CartesianCoordinates *pc,
	      CartesianCoordinates *pu, CartesianCoordinates *pv,
	      const CentralBody *pcb)
{
  SphericalCoordinates p;
  CartesianCoordinates w;

/* Find the geocentric equatorial spherical coordinates of the specified */
/* longitude at the current time.  Then rotate all vectors by minus this */
/* angle about the z-axis.                                               */

  lat_lon_to_spherical (0.0, (double) Longitude_Center_Line, t, pcb, &p);

  rotate_z(pc, -p.theta, &w);
  pc->x = w.x;
  pc->y = w.y;
  pc->z = w.z;

  rotate_z(pu, -p.theta, &w);
  pu->x = w.x;
  pu->y = w.y;
  pu->z = w.z;

  rotate_z(pv, -p.theta, &w);
  pv->x = w.x;
  pv->y = w.y;
  pv->z = w.z;
}

/*
 * Given the normalized distance to the spacecraft, compute the radius of
 * a circular footprint.
 */
static double
footprint_circle_radius( double L, int fp_angle_type, double footprint_angle )
{
  double sa = sin(footprint_angle);
  double ca = cos(footprint_angle);

  if (fp_angle_type == MASK_ELEVATION) {

    if (debug) fprintf(stderr, "Mask elevation = %f radians.\n", footprint_angle);

    return( ca*(sqrt(1-ca*ca/L/L)-sa/L) );
  } else {

    if (debug) fprintf(stderr, "Sat cone angle = %f radians.\n", footprint_angle);

    if ( L >= 1/sa )
      return( sqrt(L*L-1)/L );
    else
      return( sa*(L*ca-sqrt(1-L*L*sa*sa)) );
  }
}

/*
 * Increment the intensity image, circular footprint.
 */
static void
intensity_circle_footprint(CartesianCoordinates *pc, CartesianCoordinates *pu,
			   CartesianCoordinates *pv, int current_projection,
			   grid *g)
{
  double t, dt, t_final;
  SphericalCoordinates point;
  int j, row, prev_row, left[2], right[2], left_edge[2], right_edge[2];

  /* compute phi increment so that we fill in each row of grid */
  switch (current_projection) {
  case UNPROJECTED:
  case SINUSOIDAL:
    dt = PI/(g->height-1)/norm(pu);
    break;
  case SPHERICAL:
  case CYLINDRICAL:
  default:
    dt = 2.0/(g->height-1)/norm(pu);
  }

  t_final = PI;

  circle_point(&point, pc, pu, pv, 0.0);
  intensity_index(left, &point, current_projection, g);
  increment_intensity(left[1], left[0], g);
  row=left[1];
  prev_row=row;

  for (t=dt; t<t_final; t+=dt) {
    circle_point(&point, pc, pu, pv, t);
    intensity_index(left, &point, current_projection, g);
    row = left[1];
    if (row != prev_row) {
      if ( abs(row-prev_row) > 1 ) {
	fprintf(stderr, "Intensity: row skipped. Prev=%i,current=%i.\n",
		prev_row, row);
      }
      circle_point(&point, pc, pu, pv, TWOPI-t);
      intensity_index(right, &point, current_projection, g);

      if (left[0] <= right[0]) {
	for (j = left[0]; j <= right[0]; j++) {
	  increment_intensity(row,j,g);
	}
      } else {
	intensity_edges(left_edge, right_edge, &point, current_projection, g);
	for(j = left_edge[0]; j <= right[0]; j++) {
	  increment_intensity(row,j,g);