/* First find zenith location on Earth of each object. */

  for (i = 1; i <= tot; i++) {
    lonz[i] = DTOR(planet[i]); latz[i] = DTOR(planetalt[i]);
    EclToEqu(&lonz[i], &latz[i]);
  }

  /* Then, convert this to local horizon altitude and azimuth. */

  for (i = 1; i <= tot; i++) if (i != _MC) {
    lonz[i] = DTOR(Mod(RTOD(lonz[_MC]-lonz[i]+lon)));
    lonz[i] = DTOR(Mod(RTOD(lonz[i]-lon+PI/2.0)));
    EquToLocal(&lonz[i], &latz[i], PI/2.0-lat);
    azi[i] = DEGREES-RTOD(lonz[i]); alt[i] = RTOD(latz[i]);
  }

  /* Now, actually print the location of each object. */

  fprintf(S,
    "Body Altitude Azimuth  Azi. Vector   %s Vector    Moon Vector\n\n",
    centerplanet ? " Sun" : " Earth");
  for (k = 1; k <= tot; k++) {
    i = k <= BASE ? k : BASE+starname[k-BASE];
    if (ignore[i] || !IsThing(i))
      continue;
    AnsiColor(objectansi[i]);
    fprintf(S, "%-4.4s: ", objectname[i]);
    PrintAltitude(alt[i]);

    /* Determine directional vector based on azimuth. */

    j = (int) (FRACT(azi[i])*60.0);
    fprintf(S, " %3d%c%02d'", (int) azi[i], DEGR1, j);
    sx = cos(DTOR(azi[i])); sy = sin(DTOR(azi[i]));
    if (dabs(sx) < dabs(sy)) {
      vx = dabs(sx / sy); vy = 1.0;
    } else {
      vy = dabs(sy / sx); vx = 1.0;
    }
    fprintf(S, " (%.2f%c %.2f%c)",
      vy, sy < 0.0 ? 's' : 'n', vx, sx > 0.0 ? 'e' : 'w');

    /* Determine distance vector of current object from Sun and Moon. */

    vx = azi[1]-azi[i]; vy = azi[2]-azi[i];
    fprintf(S, " [%6.1f%6.1f] [%6.1f%6.1f]",
      dabs(vx) < DEGHALF ? vx : Sgn(vx)*(DEGREES-dabs(vx)), alt[1]-alt[i],
      dabs(vy) < DEGHALF ? vy : Sgn(vy)*(DEGREES-dabs(vy)), alt[2]-alt[i]);
    if (i >= U_LO) {
      if (i <= U_HI)
        fprintf(S, "  Uranian #%d", i-U_LO+1);
      else
        fprintf(S, "  Star #%2d", i-S_LO+1);
    }
    printl();
  }
  AnsiColor(DEFAULT);
}


/* Display x,y,z locations of each body (in AU) with respect to the Sun */
/* (or whatever the specified center planet is), as in the -S switch.   */
/* These values were already determined when calculating the planet     */
/* positions themselves, so this procedure is basically just a loop.    */

void ChartSpace()
{
  real x, y, z;
  int i;

  fprintf(S, "Body     Angle    X axis    Y axis    Z axis    Length\n");
  for (i = 0; i <= BASE; i++) {
    if (ignore[i] || i == 2 || !IsObject(i))
      continue;
    AnsiColor(objectansi[i]);
    fprintf(S, "%c%c%c%c: ", OBJNAM(i),
      objectname[i][3] ? objectname[i][3] : ' ');
    x = spacex[i]; y = spacey[i]; z = spacez[i];
    fprintf(S, "[%7.2f] [%7.2f] [%7.3f] [%7.3f] [%7.3f]",
      planet[i], x, y, z, sqrt(x*x+y*y+z*z));
    if (i >= U_LO) {
      if (i <= U_HI)
        fprintf(S, "  Uranian #%d", i-U_LO+1);
      else
        fprintf(S, "  Star #%2d", i-S_LO+1);
    }
    printl();
  }
  AnsiColor(DEFAULT);
}


/* Print the locations of the astro-graph lines on the Earth as specified */
/* with the -L switch. This includes Midheaven and Nadir lines, zenith    */
/* positions, and locations of Ascendant and Descendant lines.            */

void ChartAstroGraph()
{
  crosstruct PTR c;
  real planet1[TOTAL+1], planet2[TOTAL+1], mc[TOTAL+1], ic[TOTAL+1],
    as[TOTAL+1], ds[TOTAL+1], as1[TOTAL+1], ds1[TOTAL+1],
    lo = Lon, longm, w, x, y, z, ad, oa, am, od, dm;
  int occurcount = 0, tot = total, i, j, k, l, m, n;

  if (exdisplay & DASHL0)
    {
    Allocate(c, sizeof(crosstruct), crosstruct PTR);
    if (c == NULL
#ifdef PC
      /* For PC's the array better not cross a segment boundary. */
      || HIWORD(LOWORD(c) + sizeof(crosstruct)) > 0
#endif
      ) {
      PrintError("Not enough memory for crossing table.");
      return;
    }
  }

#ifdef MATRIX
  for (i = 1; i <= total; i++) {
    planet1[i] = DTOR(planet[i]);
    planet2[i] = DTOR(planetalt[i]);      /* Calculate zenith location on */
    EclToEqu(&planet1[i], &planet2[i]);   /* Earth of each object.        */
  }

  /* Print header. */

  fprintf(S, "Object :");
  for (j = 0, i = 1; i <= total; i++)
    if (!ignore[i] && IsThing(i)) {
      AnsiColor(objectansi[i]);
      fprintf(S, " %c%c%c", OBJNAM(i));
      j++;
      if (column80 && j >= 17) {
        tot = i;
        i = total;
      }
    }
  AnsiColor(DEFAULT);
  fprintf(S, "\n------ :");
  for (i = 1; i <= tot; i++)
    if (!ignore[i] && IsThing(i))
      fprintf(S, " ###");

  /* Print the longitude locations of the Midheaven lines. */

  fprintf(S, "\nMidheav: ");
  if (lo < 0.0)
    lo += DEGREES;
  for (i = 1; i <= tot; i++)
    if (!ignore[i] && IsThing(i)) {
    AnsiColor(objectansi[i]);
    x = DTOR(MC)-planet1[i];
    if (x < 0.0)
      x += 2.0*PI;
    if (x > PI)
      x -= 2.0*PI;
    z = lo+RTOD(x);
    if (z > DEGHALF)
      z -= DEGREES;
    mc[i] = z;
    fprintf(S, "%3.0f%c", dabs(z), z < 0.0 ? 'e' : 'w');
  }
  AnsiColor(DEFAULT);

  /* The Nadir lines are just always 180 degrees away from the Midheaven. */

  fprintf(S, "\nNadir  : ");
  for (i = 1; i <= tot; i++)
    if (!ignore[i] && IsThing(i)) {
    AnsiColor(objectansi[i]);
    z = mc[i] + DEGHALF;
    if (z > DEGHALF)
      z -= DEGREES;
    ic[i] = z;
    fprintf(S, "%3.0f%c", dabs(z), z < 0.0 ? 'e' : 'w');
  }
  AnsiColor(DEFAULT);

  /* Print the Zenith latitude locations. */

  fprintf(S, "\nZenith : ");
  for (i = 1; i <= tot; i++)
    if (!ignore[i] && IsThing(i)) {
      AnsiColor(objectansi[i]);
      y = RTOD(planet2[i]);
      fprintf(S, "%3.0f%c", dabs(y), y < 0.0 ? 's' : 'n');
      as[i] = ds[i] = as1[i] = ds1[i] = LARGE;
    }
  printl2();

  /* Now print the locations of Ascendant and Descendant lines. Since these */
  /* are curvy, we loop through the latitudes, and for each object at each  */
  /* latitude, print the longitude location of the line in question.        */

  longm = DTOR(Mod(MC+lo));
  for (j = 80; j >= -80; j -= graphstep) {
    AnsiColor(DEFAULT); 
    fprintf(S, "Asc@%2d%c: ", j >= 0 ? j : -j, j < 0 ? 's' : 'n');
    for (i = 1; i <= tot; i++)
      if (!ignore[i] && IsThing(i)) {
      AnsiColor(objectansi[i]);
      ad = tan(planet2[i])*tan(DTOR(j));
      if (ad*ad > 1.0) {
        fprintf(S, " -- ");
        as1[i] = ds1[i] = ret2[i] = LARGE;
      } else {
        ad = ASIN(ad);
        oa = planet1[i]-ad;
        if (oa < 0.0)
          oa += 2.0*PI;
        am = oa-PI/2.0;
        if (am < 0.0)
          am += 2.0*PI;
        z = longm-am;
        if (z < 0.0)
          z += 2.0*PI;
        if (z > PI)
          z -= 2.0*PI;
        as1[i] = as[i];
        as[i] = z = RTOD(z);
        ret2[i] = ad;
        fprintf(S, "%3.0f%c", dabs(z), z < 0.0 ? 'e' : 'w');
      }
    }

    /* Again, the Descendant position is related to the Ascendant's,  */
    /* being a mirror image, so it can be calculated somewhat easier. */

    AnsiColor(DEFAULT);
    fprintf(S, "\nDsc@%2d%c: ", j >= 0 ? j : -j, j < 0 ? 's' : 'n');
    for (i = 1; i <= tot; i++)
      if (!ignore[i] && IsThing(i)) {
      AnsiColor(objectansi[i]);
      ad = ret2[i];
      if (ad == LARGE)
        fprintf(S, " -- ");
      else {
        od = planet1[i]+ad;
        dm = od+PI/2.0;
        z = longm-dm;
        if (z < 0.0)
          z += 2.0*PI;
        if (z > PI)
          z -= 2.0*PI;
        ds1[i] = ds[i];
        ds[i] = z = RTOD(z);
        fprintf(S, "%3.0f%c", dabs(z), z < 0.0 ? 'e' : 'w');
      }
    }
    printl();
#endif /* MATRIX */

    /* Now, if the -L0 switch is in effect, then take these line positions, */
    /* which we saved in an array above as we were printing them, and       */
    /* calculate and print the latitude crossings.                          */

    if (exdisplay & DASHL0)
      for (l = 1; l <= total; l++) if (!ignore[l] && IsThing(l))
        for (k = 1; k <= total; k++) {
          if (ignore[k] || !IsThing(k))
            continue;
          for (n = 0; n <= 1; n++) {
            x = n ? ds1[l] : as1[l];
            y = n ? ds[l] : as[l];
            for (m = 0; m <= 1; m++) {

            /* Check if Ascendant/Descendant cross Midheaven/Nadir. */

            z = m ? ic[k] : mc[k];
            if (occurcount < MAXCROSS &&
              dabs(x-y) < DEGHALF && Sgn(z-x) != Sgn(z-y)) {
              c->obj1[occurcount] = n ? -l : l;
              c->obj2[occurcount] = m ? -k : k;
              c->lat[occurcount] = (real)j+5.0*dabs(z-y)/dabs(x-y);
              c->lon[occurcount] = z;
              occurcount++;
            }

            /* Check if Ascendant/Descendant cross another Asc/Des. */

            w = m ? ds1[k] : as1[k];
            z = m ? ds[k] : as[k];
            if (occurcount < MAXCROSS && k > l &&
                dabs(x-y)+dabs(w-z) < DEGHALF && Sgn(w-x) != Sgn(z-y)) {
              c->obj1[occurcount] = n ? -l : l;
              c->obj2[occurcount] = 100+(m ? -k : k);
              c->lat[occurcount] = (real)j+5.0*
                dabs(y-z)/(dabs(x-w)+dabs(y-z));
              c->lon[occurcount] = MIN(x, y)+dabs(x-y)*
                dabs(y-z)/(dabs(x-w)+dabs(y-z));
              occurcount++;
            }
          }
        }
    }
  }
  if ((exdisplay & DASHL0) == 0)
    return;
  printl();

  /* Now, print out all the latitude crossings we found.  */
  /* First, we sort them in order of decreasing latitude. */

  for (i = 1; i < occurcount; i++) {
    j = i-1;
    while (j >= 0 && c->lat[j] < c->lat[j+1]) {
      SWAP(c->obj1[j], c->obj1[j+1]); SWAP(c->obj2[j], c->obj2[j+1]);
      SwapReal(&c->lat[j], &c->lat[j+1]); SwapReal(&c->lon[j], &c->lon[j+1]);
      j--;
    }
  }
  for (i = 1; i < occurcount; i++) {
    j = abs(c->obj1[i]);
    AnsiColor(objectansi[j]);
    fprintf(S, "%c%c%c ", OBJNAM(j));
    AnsiColor(elemansi[c->obj1[i] > 0 ? 0 : 2]);
    fprintf(S, "%s ", c->obj1[i] > 0 ? "Ascendant " : "Descendant");
    AnsiColor(WHITE);
    fprintf(S, "crosses ");
    j = abs(c->obj2[i] - (c->obj2[i] < 50 ? 0 : 100));
    AnsiColor(objectansi[j]);
    fprintf(S, "%c%c%c ", OBJNAM(j));
    AnsiColor(elemansi[c->obj2[i] < 50 ?
      (c->obj2[i] > 0 ? 1 : 3) : (c->obj2[i] > 100 ? 0 : 2)]);
    fprintf(S, "%s ", c->obj2[i] < 50 ? (c->obj2[i] > 0 ? "Midheaven " :
      "Nadir     ") : (c->obj2[i] > 100 ? "Ascendant " : "Descendant"));
    j = (int) (FRACT(dabs(c->lon[i]))*60.0);
    AnsiColor(DEFAULT);
    fprintf(S, "at %3d%c%02d'%c, ", (int) dabs(c->lon[i]), DEGR2,
      j, c->lon[i] < 0.0 ? 'E' : 'W');
    j = (int) (FRACT(dabs(c->lat[i]))*60.0);
    fprintf(S, "%2d%c%02d'%c\n", (int) dabs(c->lat[i]), DEGR2,
      j, c->lat[i] < 0.0 ? 'S' : 'N');
  }
  Deallocate(c);
  if (!occurcount) {
    AnsiColor(DEFAULT);
    fprintf(S, "No latitude crossings.\n");
  }
}

/* charts.c */