/*========================================================================= StarCalc.c -- This module contains the fixed radix arithmetic to speed up the calculation of star positions. Floating point numbers are converted to the fixed radix format by multiplying by 10000, and the result is stored in a LONG integer. The routine sacrifices space (and some accuracy) in the interests of speed (the old tradeoff). The FastSin, FastCos, and FastTan functions rely on a table of 900 precalculated sines and tangents (in TrigTab.h) that are in the fixed format. Credits for Star Chart: Robert L. Hill of the Orange County, CA. Amiga Friends User Group wrote the original version of StarChart in AmigaBasic The star data and many of the main functions of this version are derived from that program. Ray R. Larson wrote the c version 1.0 of StarChart, 'intuitionizing' and enhancing the speed and functions of the original. Copyright (c) 1986 by Ray R. Larson This program may be freely distributed and copied, but may not be sold without the permission of the author. If you modify or enhance it, please include the above credits (and please send me a copy!). Ray R. Larson 6425 Central Ave. #304 El Cerrito, CA 94530 BitNet LARSON@UCBCMSA Well rrl CServe 70446,766 =========================================================================*/ /*------------Header file for all of the standard stuff----*/ /*-------------plus definitions of global structures-------*/ #include "star.h" /*------------Header file for sin/cos/tan/cotan table -----*/ #include "TrigTab.h" SHORT FirstStar, LastStar; extern struct ParamStruct Parms; extern struct Coord coords[]; /* x and y screen coordinates for stars */ extern struct star_rec StarTable[]; extern FLOAT P1,P12,P15; /* pi based values */ /* this routine fetches the sine values for angles specified in tenths of */ /* a degree (i.e. 451 is 45.1 degrees) */ LONG FastSin(tenthdegs) LONG tenthdegs; { register LONG deg, sign; if (tenthdegs == 0L) return(0L); if (tenthdegs < 0L) { sign = -1L; deg = -tenthdegs; } else { sign = 1L; deg = tenthdegs; } deg = deg % 3600L; if (deg > 1800L) { sign = (sign == -1L)? 1L : -1L; deg = deg - 1800L; } if (deg > 900L) deg = (900L-deg)+900L; return(SinTanTab[deg].sine * sign); } /* FastCos just shifts the sine curve and calls FastSin */ LONG FastCos(tenthdegs) LONG tenthdegs; { return(FastSin(tenthdegs + 900L)); } /* FastTan looks up the tangent */ LONG FastTan(tenthdegs) LONG tenthdegs; { register LONG deg, sign; if (tenthdegs == 0L) return(0L); deg = (tenthdegs < 0L) ? -tenthdegs : tenthdegs; deg = deg % 1800L; if (deg < 900L) return(SinTanTab[deg].tangent); else return(SinTanTab[(900L - deg) + 900L].tangent * -1L); } /* convert float to fixed format */ LONG ftofix(x) FLOAT x; { return((LONG)(x * 10000.0)); } /* multiple two fixed radix numbers */ LONG fixmult(x,y) LONG x, y; { long temp1, temp2; temp1 = (x * (y % 10000))/10000; temp2 = x * (y/10000); return (temp1+temp2); } /* divide two fixed radix numbers , lose a little precision here...*/ LONG fixdiv(x,y) LONG x, y; { long temp; if (x == 0L) return (0L); if (x == y) return (10000); /* i.e. 1 */ /* shift x over by 2 digits */ temp = x * 100; return (temp/y * 100); } /*************************************************************************** * FastCalc - Calculate the positions of the visible stars given the * current parms and plot them on the screen. Save positions in coords. * This version uses fixed radix arithmetic ***************************************************************************/ FCalc1() { LONG xDeg,yDeg,Xpos,Ypos, LatDeg, LatPos, LatCOS, LatSIN, F, LST; SHORT i; LONG yeardif, RA, DC; BOOL visible; struct star_rec *ST; struct Coord *c; /* set some (variable) constants */ yeardif = Parms.Year - 1950; LST = ftofix(Parms.Lst); LatDeg = (LONG)(Parms.Latitude * 10.0); LatPos = 1600000 - fixmult(17800, ftofix(Parms.Latitude)); LatCOS = FastCos(LatDeg); LatSIN = -FastSin(LatDeg); F = 891268; /* ie 140/(PI/2) * 10000 */ FirstStar = 0; /* for a little extra speed we will use pointers rather than array indexes */ ST = &StarTable[1]; c = &coords[1]; for (i = 1; i <= NumStars; i++) { if (yeardif) { RA = ftofix(ST->Ra); DC = ftofix(ST->Dc); /* convert right Asc. and declination of star to tenth Degrees */ xDeg = (RA * 15)/1000; yDeg = DC/1000; /* correct for year - uses mathffp routines*/ Xpos = RA + (30700 + fixmult(fixmult(13400, FastSin(xDeg)), tan(xDeg))) * yeardif / 3600; Ypos = DC + fixmult(200000, FastCos(yDeg)) * yeardif / 3600; } else { /* year is 1950 - basis for the star table */ Xpos = RA; Ypos = DC; } /* modify the Xpos and Ypos for the date, time, location and horizon */ if(Parms.Horizon == 0) visible = FastNorth(&Xpos,&Ypos,LatDeg,LatPos,LatCOS,LatSIN,F,LST); else visible = FastSouth(&Xpos,&Ypos,LatDeg,LatCOS,LatSIN,LST); /* if the star is visible from the current location - save coords and */ /* plot it on the display (not any more - actual plot is done in redisp */ if(visible) { /* set the display coordinates */ c->x = CHARTLEFT + (((2 * Xpos) + 5000)/10000); c->y = (CHARTTOP+1L) + ((Ypos + 5000)/10000); /* set the first and last coords filled in */ if (FirstStar) LastStar = i; else FirstStar = i; } else { c->x = 0L; c->y = 0L; } /* increment the pointers */ ST++; c++; } /* end of for loop */ } /* end of DisplayStars */ /*************************************************************************** * FastNorth - calculate star position for Northern Sky * returns 1 if star is visible, zero if not visible from current location ***************************************************************************/ BOOL FastNorth (xposit,yposit,lat,LatPos,latcos,latsin,F,LST) LONG *xposit, *yposit,lat,LatPos,latcos,latsin,F,LST; { LONG DecDeg,HourDeg, Temp1, TempX, TempY; register LONG xpos, ypos; xpos = *xposit; ypos = *yposit; if (ypos < 0) return(0); /* convert declination, latitude, and local siderial time to Degrees */ DecDeg = ypos/1000; HourDeg = (LST - xpos)*15 / 1000; /* eliminate stars below the horizon line */ if(( fixmult(latcos,fixmult(FastCos(DecDeg), FastCos(HourDeg)))) <( fixmult(latsin,FastSin(DecDeg)))) return(0); xpos = xpos - LST; /* Map to scale for the display */ xpos = (xpos * 15) - 15708; Temp1 = fixmult(F, (15708 - abs(ypos))); TempX = fixmult(Temp1,FastCos((xpos/1000))) + 1400000; TempY = fixmult(Temp1, FastSin((xpos/1000))) + LatPos; xpos = TempX; ypos = TempY; if((xpos > 2790000) || (xpos < 0)) return(0); if((ypos > 1590000) || (ypos < 0)) return(0); /* if we got to here it must be visible */ *xposit = xpos; *yposit = ypos; return(1); } /* end of FastNorth */ /*************************************************************************** * FastSouth - calculate star position for Southern Sky * returns 1 if star is visible, zero if not visible from current location ***************************************************************************/ BOOL FastSouth (xposit,yposit,latdeg,latcos,latsin,LST) LONG *xposit, *yposit,latdeg,latcos,latsin,LST; { LONG DecDeg,LatDeg,HourDeg, LatPos, Temp1; register LONG xpos, ypos; xpos = *xposit; ypos = *yposit; if (ypos > (latdeg*1000)) return(0); if (ypos < ((latdeg*1000) - 900000)) return(0); /* convert declination, latitude, and local siderial time to Degrees */ DecDeg = ypos/1000; HourDeg = (LST - xpos)*15 / 1000; /* eliminate stars below the horizon line */ if(( fixmult(latcos,fixmult(FastCos(DecDeg), FastCos(HourDeg)))) <( fixmult(latsin,FastSin(DecDeg)))) return(0); xpos = LST - xpos; if(xpos < -120000) xpos += 240000; if(xpos > 120000) xpos -= 240000; xpos = 1400000 + fixmult(233300 ,xpos); if((xpos > 2790000) || (xpos < 0)) return(0); ypos = fixmult(17800, (latdeg*1000) - ypos); if((ypos > 1590000) || (ypos < 0)) return(0); /* if we got to here it must be visible */ *xposit = xpos; *yposit = ypos; return(1); } /* end of FastSouth */ FLOAT FSin(x) FLOAT x; { LONG tenthdegs; tenthdegs = (LONG) (x * 57.2958)* 10.0; /* convert to degrees */ return(((FLOAT)FastSin(tenthdegs)/ 10000.0)); } FLOAT FCos(x) FLOAT x; { LONG tenthdegs; tenthdegs = (LONG) (x * 57.2958)* 10.0; /* convert to degrees */ return(((FLOAT)FastCos(tenthdegs)/ 10000.0)); } FLOAT FTan(x) FLOAT x; { LONG tenthdegs; tenthdegs = (LONG) (x * 57.2958)* 10.0; /* convert to degrees */ return(((FLOAT)FastTan(tenthdegs)/ 10000.0)); } /*************************************************************************** * CalcStars - Calculate the positions of the visible stars given the * current parms and plot them on the screen. Save positions in coords. ***************************************************************************/ FastCalc() { FLOAT xRad,yRad,Xpos,Ypos, LatRad, LatPos, LatCOS, LatSIN, F; SHORT i; FLOAT yeardif, RA, DC; BOOL visible; struct star_rec *ST; struct Coord *c; /* set some (variable) constants */ yeardif = (FLOAT)(Parms.Year - 1950); LatRad = Parms.Latitude * P1; LatPos = 160.0 - (1.78 * Parms.Latitude); LatCOS = FCos(LatRad); LatSIN = -FSin(LatRad); F = 140.0 / 1.5708; /* for a little extra speed we will use pointers rather than array indexes */ ST = &StarTable[1]; c = &coords[1]; for (i = 1; i <= NumStars; i++) { if (((Parms.Horizon == 0)&&(ST->Dc < 0.0)) || /* north view */ ((Parms.Horizon == 1)&& ((ST->Dc > (FLOAT)Parms.Latitude)|| (ST->Dc < (FLOAT)(Parms.Latitude - 90))) /* south view */ )) { /* don't bother calculating if it will never be visible given the */ /* current position and orientation. */ c->x = 0L; c->y = 0L; ST++; c++; continue; } if (yeardif) { RA = ST->Ra; DC = ST->Dc; /* convert right Asc. and declination of star to radians */ xRad = RA * P12; yRad = DC * P1; /* correct for year - uses mathffp routines*/ Xpos = RA + (3.07 + 1.34 * FSin(xRad) * FTan(xRad)) * yeardif / 3600.0; Ypos = DC + 20.0 * FCos(yRad) * yeardif / 3600.0; } else { /* year is 1950 - basis for the star table */ Xpos = RA; Ypos = RA; } /* modify the Xpos and Ypos for the date, time, location and horizon */ if(Parms.Horizon == 0) visible = NorthVF(&Xpos,&Ypos,LatRad,LatPos,LatCOS,LatSIN,F); else visible = SouthVF(&Xpos,&Ypos,LatCOS,LatSIN); /* if the star is visible from the current location - save coords and */ /* plot it on the display (not any more - actual plot is done in redisp */ if(visible) { /* set the display coordinates */ c->x = CHARTLEFT + (LONG)((2.0 * Xpos) + 0.5); c->y = (CHARTTOP+1L) + (LONG)(Ypos + 0.5); } else { c->x = 0L; c->y = 0L; } /* increment the pointers */ ST++; c++; } /* end of for loop */ } /* end of DisplayStars */ /*************************************************************************** * NorthView - calculate star position for Northern Sky * returns 1 if star is visible, zero if not visible from current location ***************************************************************************/ BOOL NorthVF(xposit,yposit,lat,LatPos,latcos,latsin,F) FLOAT *xposit, *yposit,lat,LatPos,latcos,latsin,F; { FLOAT DecRad,HourRad, Temp1, TempX, TempY; register FLOAT xpos, ypos; xpos = *xposit; ypos = *yposit; if (ypos < 0.0) return(0); /* convert declination, latitude, and local siderial time to radians */ DecRad = ypos * P1; HourRad = P12 * (Parms.Lst - xpos); /* eliminate stars below the horizon line */ if(( latcos * FCos(DecRad) * FCos(HourRad)) <( latsin * FSin(DecRad))) return(0); xpos = xpos - Parms.Lst; ypos = ypos * P1; /* Map to scale for the display */ xpos = (xpos * P15) - 1.5708; Temp1 = F * (1.5708 - abs(ypos)); TempX = Temp1 * FCos(xpos) + 140.0; TempY = Temp1 * FSin(xpos) + LatPos; xpos = TempX; ypos = TempY; if((xpos > 279.0) || (xpos < 0.0)) return(0); if((ypos > 159.0) || (ypos < 0.0)) return(0); /* if we got to here it must be visible */ *xposit = xpos; *yposit = ypos; return(1); } /* end of NorthView */ /*************************************************************************** * SouthView - calculate star position for Southern Sky * returns 1 if star is visible, zero if not visible from current location ***************************************************************************/ BOOL SouthVF (xposit,yposit,latcos,latsin) FLOAT *xposit, *yposit,latcos,latsin; { FLOAT DecRad,LatRad,HourRad, LatPos, Temp1; register FLOAT xpos, ypos; xpos = *xposit; ypos = *yposit; if (ypos > (FLOAT)Parms.Latitude) return(0); if (ypos < (FLOAT)(Parms.Latitude - 90)) return(0); /* convert declination, latitude, and local siderial time to radians */ DecRad = ypos * P1; HourRad = P12 * (Parms.Lst - xpos); /* eliminate stars below the horizon line */ if(( latcos * FCos(DecRad) * FCos(HourRad)) <( latsin * FSin(DecRad))) return(0); xpos = Parms.Lst - xpos; if(xpos < -12.0) xpos += 24.0; if(xpos > 12.0) xpos -= 24.0; xpos = 140.0 + (23.33 * xpos); if((xpos > 279) || (xpos < 0)) return(0); ypos = 1.78 * (Parms.Latitude - ypos); if((ypos > 159.0) || (ypos < 0.0)) return(0); /* if we got to here it must be visible */ *xposit = xpos; *yposit = ypos; return(1); } /* end of SouthView */