/* Now display each aspect line. */

  for (i = 0; i < occurcount; i++) {
    fprintf(S, "%3d: ", i+1);
    j = source[i]; k = aspect[i]; l = dest[i];
    PrintAspect(
      j, ZTOS(planet[j]), (int)Sgn(ret[j]), k,
      l, ZTOS(planet[l]), (int)Sgn(ret[l]), 'D');
    m = grid->v[j][l];
    AnsiColor(m < 0 ? WHITE : LTGRAY);
    fprintf(S, "- %s%2d%c%02d'", m < 0 ? "app" : "sep",
      abs(m)/60, DEGR1, abs(m)%60);
    AnsiColor(DKGREEN);
    fprintf(S, " - power:%6.2f", power[i]);
    PrintInDay(j, k, l);
  }
}


/* Search through a month, and print out the times of exact transits where   */
/* planets in the time frame make aspect to the planets in some other chart, */
/* as specified with the -t switch. To do this, we cast charts for the start */
/* and end of the month, or part of the month, and do an equation check for  */
/* aspects to the other base chart during the interval.                      */

void DisplayTransitSearch(prog)
int prog;
{
  real planet3[TOTAL+1], house3[SIGNS+1], ret3[TOTAL+1];
  uint time[MAXINDAY];
  int source[MAXINDAY], aspect[MAXINDAY], dest[MAXINDAY], sign[MAXINDAY],
    isret[MAXINDAY], M1, M2, occurcount, div, i, j, k, s1, s2, s3;
  real divsiz, daysiz, d, e1, e2, f1, f2;

  for (i = 1; i <= SIGNS; i++)
    house3[i] = house[i];
  for (i = 1; i <= total; i++) {
    planet3[i] = planet[i];
    ret3[i] = ret[i];
  }
  if (Mon2 == 0) {      /* Searching month number zero means to search */
    M1 = 1; M2 = 12;    /* the whole year instead, month by month.     */
  } else
    M1 = M2 = Mon2;

  /* Start searching the month or months in question for any transits. */

  for (Mon2 = M1; Mon2 <= M2; Mon2++) {
    daysiz = (real) DayInMonth(Mon2, Yea2)*24.0*60.0;
    divsiz = daysiz/ (real) divisions;

    /* Cast chart for beginning of month and store it for future use. */

    SetCore(Mon2, 1, Yea2, 0.0, Zon2, Lon2, Lat2);
    if (progress = prog) {
      Jdp = (real)MdyToJulian(MM, DD, YY);
      SetCore(Mon, Day, Yea, Tim, Zon, Lon, Lat);
    }
    CastChart(TRUE);
    for (i = 1; i <= SIGNS; i++)
      house2[i] = house[i];
    for (i = 1; i <= OBJECTS; i++) {
      planet2[i] = planet[i];
      ret2[i] = ret[i];
    }

    /* Divide our month into segments and then search each segment in turn. */

    for (div = 1; div <= divisions; div++) {
      occurcount = 0;

      /* Cast the chart for the ending time of the present segment, and */
      /* copy the start time chart from the previous end time chart.    */

      d = 1.0 + (daysiz/24.0/60.0)*(real)div/(real)divisions;
      SetCore(Mon2, (int)d, Yea2,
        DegToDec(FRACT(d)*24.0), Zon2, Lon2, Lat2);
      if (prog) {
        Jdp = (real)MdyToJulian(MM, DD, YY);
        SetCore(Mon, Day, Yea, Tim, Zon, Lon, Lat);
      }
      CastChart(TRUE);
      for (i = 1; i <= SIGNS; i++) {
        house1[i] = house2[i]; house2[i] = house[i];
      }
      for (i = 1; i <= OBJECTS; i++) {
        planet1[i] = planet2[i]; ret1[i] = ret2[i];
        planet2[i] = planet[i];  ret2[i] = ret[i];
      }

      /* Now search through the present segment for any transits. Note that */
      /* stars can be transited, but they can't make transits themselves.   */

      for (i = 1; i <= total; i++) if (!ignore[i])
        for (j = 1; j <= BASE; j++) if (!ignore2[j] && (prog || IsThing(j)))

          /* Between each pair of planets, check if they make any aspects. */

          for (k = 1; k <= aspects; k++) {
            d = planet3[i]; e1 = planet1[j]; e2 = planet2[j];
            if (MinDistance(e1, Mod(d-aspectangle[k])) <
                MinDistance(e2, Mod(d+aspectangle[k]))) {
              e1 = Mod(e1+aspectangle[k]);
              e2 = Mod(e2+aspectangle[k]);
            } else {
              e1 = Mod(e1-aspectangle[k]);
              e2 = Mod(e2-aspectangle[k]);
            }

            /* Check to see if the present aspect actually occurs during the */
            /* segment, making sure we check any Aries point crossings.      */

            f1 = e1-d;
            if (dabs(f1) > DEGHALF)
              f1 -= Sgn(f1)*DEGREES;
            f2 = e2-d;
            if (dabs(f2) > DEGHALF)
              f2 -= Sgn(f2)*DEGREES;
            if (MinDistance(d, Midpoint(e1, e2)) < DEGQUAD &&
              Sgn(f1) != Sgn(f2) && occurcount < MAXINDAY) {

              /* Ok, we have found a transit. Now determine the time */
              /* and save this transit in our list to be printed.    */

              source[occurcount] = j;
              aspect[occurcount] = k;
              dest[occurcount] = i;
              time[occurcount] = (int) (dabs(f1)/(dabs(f1)+dabs(f2))*divsiz) +
                (int) ((real) (div-1)*divsiz);
              sign[occurcount] = (int) (Mod(
                MinDistance(planet1[j], Mod(d-aspectangle[k])) <
                MinDistance(planet2[j], Mod(d+aspectangle[k])) ?
                d-aspectangle[k] : d+aspectangle[k])/30.0)+1;
              isret[occurcount] = (int)Sgn(ret1[j]) + (int)Sgn(ret2[j]);
              occurcount++;
            }
          }

      /* After all transits located, sort them by time at which they occur. */

      for (i = 1; i < occurcount; i++) {
        j = i-1;
        while (j >= 0 && time[j] > time[j+1]) {
          SWAP(source[j], source[j+1]);
          SWAP(aspect[j], aspect[j+1]);
          SWAP(dest[j], dest[j+1]);
          SWAP(time[j], time[j+1]);
          SWAP(sign[j], sign[j+1]);
          SWAP(isret[j], isret[j+1]);
          j--;
        }
      }

      /* Now loop through list and display all the transits. */

      for (i = 0; i < occurcount; i++) {
        k = smartcusp && IsCusp(dest[i]);
        if (k && aspect[i] > _OPP)
          continue;
        else
          k = k && aspect[i] == _OPP;
        s1 = time[i]/24/60;
        s3 = time[i]-s1*24*60;
        s2 = s3/60;
        s3 = s3-s2*60;
        SetSave(Mon2, s1+1, Yea2,
          DegToDec((real)(time[i]-s1*24*60) / 60.0), Zon2, Lon2, Lat2);
        fprintf(S, "%s %s ",
          CharDate(Mon2, s1+1, Yea2, FALSE), CharTime(s2, s3));
        PrintAspect(source[i], sign[i], isret[i], aspect[i],
          dest[i], ZTOS(planet3[dest[i]]), (int)Sgn(ret3[dest[i]]),
          prog ? 'u' : 't');

        /* Check for a Solar, Lunar, or any other return. */

        if (aspect[i] == _CON && source[i] == dest[i]) {
          AnsiColor(WHITE);
          fprintf(S, " (%s Return)", source[i] == _SUN ? "Solar" :
            (source[i] == _MOO ? "Lunar" : objectname[source[i]]));
        }
        printl();
#ifdef INTERPRET
        if (interpret)
          InterpretTransit(source[i], aspect[i], dest[i]);
#endif
        AnsiColor(DEFAULT);
      }
    }
  }
}


/* Given an arbitrary day, determine what aspects are made between this */
/* transiting chart and the given natal chart, as specified with the -T */
/* switch, and display the transits in order sorted by influence.       */

void DisplayTransitInfluence(prog)
int prog;
{
  int source[MAXINDAY], aspect[MAXINDAY], dest[MAXINDAY];
  real power[MAXINDAY];
  byte ignore3[TOTAL+1];
  int occurcount = 0, i, j, k, l, m;

  /* Cast the natal and transiting charts as with a relationship chart. */

  for (i = 1; i <= SIGNS; i++)
    house1[i] = house[i];
  for (i = 1; i <= total; i++) {
    planet1[i] = planet[i];
    ret1[i] = ret[i];
  }
  SetCore(Mon2, Day2, Yea2, Tim2, Zon2, Lon2, Lat2);
  if (progress = prog) {
    Jdp = (real)MdyToJulian(MM, DD, YY);
    MM = Mon; DD = Day; YY = Yea;
  }
  CastChart(TRUE);
  for (i = 1; i <= SIGNS; i++)
    house2[i] = house[i];
  for (i = 1; i <= total; i++) {
    planet2[i] = planet[i];
    ret2[i] = ret[i];
  }

  /* Do a relationship aspect grid to get the transits. We have to make and */
  /* restore three changes to get it right for this chart. (1) We make the  */
  /* natal planets have zero velocity so applying vs. separating is only a  */
  /* function of the transiter. (2) We force applying vs. separating orbs   */
  /* regardless if -ga or -ma is in effect or not. (3) Finally we tweak the */
  /* main restrictions to allow for transiting objects not restricted.      */

  for (i = 1; i <= total; i++) {
    ret[i] = ret1[i];
    ret1[i] = 0.0;
    ignore3[i] = ignore[i];
    ignore[i] = ignore[i] && ignore2[i];
  }
  i = exdisplay;
  exdisplay |= DASHga;
  CreateGridRelation(FALSE);
  exdisplay = i;
  for (i = 1; i <= total; i++) {
    ret1[i] = ret[i];
    ignore[i] = ignore3[i];
  }

  /* Loop through the grid, and build up a list of the valid transits. */

  for (i = 1; i <= BASE; i++) {
    if (ignore2[i] || !IsThing(i))
      continue;
    for (j = 1; j <= total; j++) {
      if (ignore[j] || (k = grid->n[i][j]) == 0 || occurcount >= MAXINDAY)
        continue;
      if (smartcusp && k > _OPP && IsCusp(j))
        continue;
      source[occurcount] = i; aspect[occurcount] = k; dest[occurcount] = j;
      l = grid->v[i][j];
      power[occurcount] = transitinf[i]*
        ((j <= BASE ? objectinf[j] : 2.0)/4.0)*aspectinf[k]*
        (1.0-(real)abs(l)/60.0/Orb(i, j, k));
      occurcount++;
    }
  }

  /* After all transits located, sort them by their total power. */

  for (i = 1; i < occurcount; i++) {
    j = i-1;
    while (j >= 0 && power[j] < power[j+1]) {
      SWAP(source[j], source[j+1]);
      SWAP(aspect[j], aspect[j+1]);
      SWAP(dest[j], dest[j+1]);
      SwapReal(&power[j], &power[j+1]);
      j--;
    }
  }

  /* Now loop through list and display each transit in effect at the time. */

  for (i = 0; i < occurcount; i++) {
    k = aspect[i];
    l = source[i];
    fprintf(S, "%3d: ", i+1);
    j = ZTOS(planet2[l]);
    PrintAspect(l, j, (int)Sgn(ret2[l]), k,
      dest[i], ZTOS(planet1[dest[i]]), (int)Sgn(ret1[dest[i]]),
      prog ? 'U' : 'T');
    m = grid->v[l][dest[i]];
    AnsiColor(m < 0 ? WHITE : LTGRAY);
    fprintf(S, "- %s%2d%c%02d'", m < 0 ? "app" : "sep",
      abs(m)/60, DEGR1, abs(m)%60);
    AnsiColor(DKGREEN);
    fprintf(S, " - power:%6.2f", power[i]);
    if (k == _CON && l == dest[i]) {       /* Print a small "R" for returns. */
      AnsiColor(WHITE);
      fprintf(S, " R");
    }
    printl();
#ifdef INTERPRET
    if (interpret)
      InterpretTransit(l, k, dest[i]);
#endif
    AnsiColor(DEFAULT);
  }
}


/* Given the zodiac location of a planet in the sky and its declination,   */
/* and a location on the Earth, compute the azimuth and altitude of where  */
/* on the local horizon sky the planet would appear to one at the given    */
/* location. A reference MC position at Greenwich is also needed for this. */

void EclToHorizon(azi, alt, planet, planetalt, lon, lat, mc)
real *azi, *alt, planet, planetalt, lon, lat, mc;
{
  real lonz, latz;

  lonz = DTOR(planet); latz = DTOR(planetalt);
  EclToEqu(&lonz, &latz);
  lonz = DTOR(Mod(RTOD(mc-lonz+lon)));
  lonz = DTOR(Mod(RTOD(lonz-lon+PI/2.0)));
  EquToLocal(&lonz, &latz, PI/2.0-lat);
  *azi = DEGREES-RTOD(lonz); *alt = RTOD(latz);
}


/* Display a list of planetary rising times relative to the local horizon */
/* for the day indicated in the chart information, as specified with the  */
/* -Zd switch. For the day, the time each planet rises (transits horizon  */
/* in East half of sky), sets (transits horizon in West), reaches its     */
/* zenith point (transits meridian in South half of sky), and nadirs      */
/* transits meridian in North), is displayed.                             */

void DisplayInDayHorizon()
{
  int source[MAXINDAY], type[MAXINDAY], time[MAXINDAY], sign[MAXINDAY],
    isret[MAXINDAY], occurcount, division, div, s1, s2, i, j;
  real planetalt1[TOTAL+1], planetalt2[TOTAL+1], azialt[MAXINDAY],
    divsiz, lon, lat, mc1, mc2, azi1, alt1, azi2, alt2, d, k;

  lon = DTOR(Mod(Lon)); lat = DTOR(Lat);
  division = divisions * 4;
  divsiz = 24.0/ (real) division*60.0;
  occurcount = 0;

  SetCore(Mon, Day, Yea, 0.0, Zon, Lon, Lat);
  CastChart(TRUE);
  mc2 = DTOR(planet[_MC]); k = DTOR(planetalt[_MC]);
  EclToEqu(&mc2, &k);
  for (i = 1; i <= SIGNS; i++) {
    house2[i] = house[i];
    inhouse2[i] = inhouse[i];
  }
  for (i = 1; i <= total; i++) {
    planet2[i] = planet[i];
    planetalt2[i] = planetalt[i];
    ret2[i] = ret[i];
  }

  /* Loop thorough the day, dividing it into a certain number of segments. */
  /* For each segment we get the planet positions at its endpoints.        */

  for (div = 1; div <= division; div++) {
    SetCore(Mon, Day, Yea,
      DegToDec(24.0*(real)div/(real)division), Zon, Lon, Lat);
    CastChart(TRUE);
    mc1 = mc2;
    mc2 = DTOR(planet[_MC]); k = DTOR(planetalt[_MC]);
    EclToEqu(&mc2, &k);
    for (i = 1; i <= SIGNS; i++) {
      house1[i] = house2[i]; inhouse1[i] = inhouse2[i];
      house2[i] = house[i];  inhouse2[i] = inhouse[i];
    }
    for (i = 1; i <= total; i++) {
      planet1[i] = planet2[i]; planet2[i] = planet[i];
      planetalt1[i] = planetalt2[i]; planetalt2[i] = planetalt[i];
      ret1[i] = ret2[i]; ret2[i] = ret[i];
    }

    /* For our segment, check to see if each planet during it rises, sets, */
    /* reaches its zenith, or reaches its nadir.                           */

    for (i = 1; i <= total; i++) if (!ignore[i] && IsThing(i)) {
      EclToHorizon(&azi1, &alt1, planet1[i], planetalt1[i], lon, lat, mc1);
      EclToHorizon(&azi2, &alt2, planet2[i], planetalt2[i], lon, lat, mc2);
      j = 0;
      /* Check for transits to the horizon. */
      if ((alt1 > 0.0) != (alt2 > 0.0)) {
        d = dabs(alt1)/(dabs(alt1)+dabs(alt2));
        k = Mod(azi1 + d*MinDifference(azi1, azi2));
        j = 1 + 2*(MinDistance(k, DEGHALF) < DEGQUAD);
      /* Check for transits to the meridian. */
      } else if (Sgn(MinDifference(azi1, DEGQUAD)) !=
        Sgn(MinDifference(azi2, DEGQUAD))) {
        j = 2 + 2*(MinDistance(azi1, DEGQUAD) < DEGQUAD);
        d = dabs(azi1 - (j > 2 ? DEGQUAD : 270.0))/MinDistance(azi1, azi2);
        k = alt1 + d*(alt2-alt1);
      }
      if (j && occurcount < MAXINDAY) {
        source[occurcount] = i;
        type[occurcount] = j;
        time[occurcount] = (int)(24.0*((real)(div-1)+d)/(real)division*60.0);
        sign[occurcount] = (int)Mod(planet1[i] +
          d*MinDifference(planet1[i], planet2[i]))/30 + 1;
        isret[occurcount] = (int)Sgn(ret1[i]) + (int)Sgn(ret2[i]);
        azialt[occurcount] = k;
        SetSave(Mon, Day, Yea, DegToDec((real)time[occurcount] / 60.0),
          Zon, Lon, Lat);
        occurcount++;
      }
    }
  }

  /* Sort each event in order of time when it happens during the day. */

  for (i = 1; i < occurcount; i++) {
    j = i-1;
    while (j >= 0 && time[j] > time[j+1]) {
      SWAP(source[j], source[j+1]);
      SWAP(type[j], type[j+1]);
      SWAP(time[j], time[j+1]);
      SWAP(sign[j], sign[j+1]);
      SWAP(isret[j], isret[j+1]);