home *** CD-ROM | disk | FTP | other *** search
/ The Starbase One Astronomy & Space Collection / STARBASE_ONE.ISO / a96 / disk07 / altaz.exe / CONVERT.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1995-08-09  |  10.0 KB  |  453 lines
  1. #include <iostream.h>
  2. #include <iomanip.h>
  3. #include <fstream.h>
  4. #include <strstrea.h>
  5. #include <math.h>
  6. #include <conio.h>
  7. #include <stdio.h>
  8. #include <string.h>
  9. #include <values.h>
  10.  
  11. #include "common.h"
  12. #include "video.hpp"
  13. #include "atimes.hpp"
  14. #include "badexit.hpp"
  15. #include "coord.hpp"
  16. #include "convert.hpp"
  17.  
  18. Convert::Convert
  19. (
  20.     const double Z1Deg,
  21.     const double Z2Deg,
  22.     const double Z3Deg,
  23.     const double LongitudeDeg,
  24.     int Tz
  25. ):
  26.     Coord(LongitudeDeg,Tz),
  27.  
  28.     Z1Deg(Z1Deg),
  29.     Z2Deg(Z2Deg),
  30.     Z3Deg(Z3Deg)
  31. {
  32.     SetMountErrorsDeg(Z1Deg,Z2Deg,Z3Deg);
  33.     InitFlag=NoneInit;
  34.     LoadOneTwo();
  35. }
  36.  
  37. void Convert::SaveOneTwo(void)
  38. {
  39.     ofstream Output(OneTwoFile);
  40.     cout.setf(ios::showpoint);
  41.     cout.setf(ios::fixed,ios::floatfield);
  42.     Output<<"One.Ra "<<setprecision(17)<<setw(17)<<One.Ra<<endl;
  43.     Output<<"One.Dec "<<setprecision(17)<<setw(17)<<One.Dec<<endl;
  44.     Output<<"One.Alt "<<setprecision(17)<<setw(17)<<One.Alt<<endl;
  45.     Output<<"One.Az "<<setprecision(17)<<setw(17)<<One.Az<<endl;
  46.     Output<<"One.SidT.Time "<<setprecision(17)<<setw(17)<<One.SidT.Time<<endl;
  47.     Output<<"Two.Ra "<<setprecision(17)<<setw(17)<<Two.Ra<<endl;
  48.     Output<<"Two.Dec "<<setprecision(17)<<setw(17)<<Two.Dec<<endl;
  49.     Output<<"Two.Alt "<<setprecision(17)<<setw(17)<<Two.Alt<<endl;
  50.     Output<<"Two.Az "<<setprecision(17)<<setw(17)<<Two.Az<<endl;
  51.     Output<<"Two.SidT.Time "<<setprecision(17)<<setw(17)<<Two.SidT.Time<<endl;
  52. }
  53.  
  54. void Convert::LoadOneTwo(void)
  55. {
  56.     char Name[NameSize];
  57.     char Value[ValueSize];
  58.     Position Temp;
  59.  
  60.     ifstream File(OneTwoFile);
  61.     if(File)
  62.     {
  63.         while(!File.eof())
  64.         {
  65.             File>>Name>>Value;
  66.             if((strncmpi(Name,"One.Ra",sizeof(Name)))==0)
  67.                 One.Ra=atof(Value);
  68.             else if((strncmpi(Name,"One.Dec",sizeof(Name)))==0)
  69.                 One.Dec=atof(Value);
  70.             else if((strncmpi(Name,"One.Alt",sizeof(Name)))==0)
  71.                 One.Alt=atof(Value);
  72.             else if((strncmpi(Name,"One.Az",sizeof(Name)))==0)
  73.                 One.Az=atof(Value);
  74.             else if((strncmpi(Name,"One.SidT.Time",sizeof(Name)))==0)
  75.                 One.SidT.Time=atof(Value);
  76.             else if((strncmpi(Name,"Two.Ra",sizeof(Name)))==0)
  77.                 Two.Ra=atof(Value);
  78.             else if((strncmpi(Name,"Two.Dec",sizeof(Name)))==0)
  79.                 Two.Dec=atof(Value);
  80.             else if((strncmpi(Name,"Two.Alt",sizeof(Name)))==0)
  81.                 Two.Alt=atof(Value);
  82.             else if((strncmpi(Name,"Two.Az",sizeof(Name)))==0)
  83.                 Two.Az=atof(Value);
  84.             else if((strncmpi(Name,"Two.SidT.Time",sizeof(Name)))==0)
  85.                 Two.SidT.Time=atof(Value);
  86.         }
  87.         Temp=Current;
  88.         Current=One;
  89.         InitMatrix(InitOne);
  90.         Current=Two;
  91.         InitMatrix(InitTwo);
  92.         Current=Temp;
  93.     }
  94.     else
  95.         cout<<endl<<OneTwoFile<<" not found";
  96. }
  97.  
  98. void Convert::SetMountErrorsDeg(const double Z1Deg,const double Z2Deg,
  99. const double Z3Deg)
  100. {
  101.     Z1=Z1Deg/RadToDeg;
  102.     Z2=Z2Deg/RadToDeg;
  103.     Z3=Z3Deg/RadToDeg;
  104.     if(Z1!=0 || Z2!=0)
  105.         Z1Z2NonZeroFlag=Yes;
  106.     else
  107.         Z1Z2NonZeroFlag=No;
  108. }
  109.  
  110. void Convert::ZeroArrays(void)
  111. {
  112.     int I,J;
  113.  
  114.     for(I=0;I<4;I++)
  115.         for(J=0;J<4;J++)
  116.             Q[I][J]=V[I][J]=R[I][J]=X[I][J]=Y[I][J]=0;
  117. }
  118.  
  119. void Convert::InitMatrix(const int Init)
  120. {
  121.     int I,J,L,M,N;
  122.     double D,B,F,H,A,E,W;
  123.  
  124.     if(Init==InitOne || (Init==InitTwo && InitFlag))
  125.     {
  126.         if(Init==InitOne)
  127.             ZeroArrays();
  128.         D=Current.Dec;
  129.         // B is CCW so HA formula backwards
  130.         B=Current.Ra-Current.SidT.Time;
  131.         X[1][Init]=cos(D)*cos(B);
  132.         X[2][Init]=cos(D)*sin(B);
  133.         X[3][Init]=sin(D);
  134.         // F is CCW
  135.         F=OneRev-Current.Az;
  136.         H=Current.Alt+Z3;
  137.         SubrA(F,H,Z1,Z2);
  138.         Y[1][Init]=Y[1][0];
  139.         Y[2][Init]=Y[2][0];
  140.         Y[3][Init]=Y[3][0];
  141.         // InitMatrix(InitOne) must already be called
  142.         if(Init==InitTwo)
  143.         {
  144.             X[1][3]=X[2][1]*X[3][2]-X[3][1]*X[2][2];
  145.             X[2][3]=X[3][1]*X[1][2]-X[1][1]*X[3][2];
  146.             X[3][3]=X[1][1]*X[2][2]-X[2][1]*X[1][2];
  147.             A=sqrt(X[1][3]*X[1][3]+X[2][3]*X[2][3]+X[3][3]*X[3][3]);
  148.             if(A==0)
  149.                 A=1/MAXDOUBLE;
  150.             for(I=1; I<=3; I++)
  151.                 X[I][3]/=A;
  152.             Y[1][3]=Y[2][1]*Y[3][2]-Y[3][1]*Y[2][2];
  153.             Y[2][3]=Y[3][1]*Y[1][2]-Y[1][1]*Y[3][2];
  154.             Y[3][3]=Y[1][1]*Y[2][2]-Y[2][1]*Y[1][2];
  155.             A=sqrt(Y[1][3]*Y[1][3]+Y[2][3]*Y[2][3]+Y[3][3]*Y[3][3]);
  156.             if(A==0)
  157.                 A=1/MAXDOUBLE;
  158.             for(I=1; I<=3; I++)
  159.                 Y[I][3]/=A;
  160.             // transform matrix
  161.             for(I=1; I<=3; I++)
  162.                 for(J=1; J<=3; J++)
  163.                     V[I][J]=X[I][J];
  164.             DeterminateSubr(W);
  165.             E=W;
  166.             if(E==0)
  167.                 E=1/MAXDOUBLE;
  168.             for(M=1; M<=3; M++)
  169.             {
  170.                 for(I=1; I<=3; I++)
  171.                     for(J=1; J<=3; J++)
  172.                         V[I][J]=X[I][J];
  173.                 for(N=1; N<=3; N++)
  174.                 {
  175.                     V[1][M]=0;
  176.                     V[2][M]=0;
  177.                     V[3][M]=0;
  178.                     V[N][M]=1;
  179.                     DeterminateSubr(W);
  180.                     Q[M][N]=W/E;
  181.                 }
  182.             }
  183.             for(I=1; I<=3; I++)
  184.                 for(J=1; J<=3; J++)
  185.                     R[I][J]=0;
  186.             for(I=1; I<=3; I++)
  187.                 for(J=1; J<=3; J++)
  188.                     for(L=1; L<=3; L++)
  189.                         R[I][J]+=(Y[I][L]*Q[L][J]);
  190.             for(M=1; M<=3; M++)
  191.             {
  192.                 for(I=1; I<=3; I++)
  193.                     for(J=1; J<=3; J++)
  194.                         V[I][J]=R[I][J];
  195.                 DeterminateSubr(W);
  196.                 E=W;
  197.                 if(E==0)
  198.                     E=1/MAXDOUBLE;
  199.                 for(N=1; N<=3; N++)
  200.                 {
  201.                     V[1][M]=0;
  202.                     V[2][M]=0;
  203.                     V[3][M]=0;
  204.                     V[N][M]=1;
  205.                     DeterminateSubr(W);
  206.                     Q[M][N]=W/E;
  207.                 }
  208.             }
  209.         }
  210.         InitFlag=Init;
  211.         if(Init==InitOne)
  212.             One=Current;
  213.         else
  214.         {
  215.             Two=Current;
  216.             SaveOneTwo();
  217.             CalcPostInitVars();
  218.         }
  219.     }
  220.     else
  221.         BadExit E("Bad Init parameter passed in Convert::InitMatrix");
  222. }
  223.  
  224. void Convert::SubrA(double F,double H,double Z1,double Z2)
  225. {
  226.     double CosF,CosH,SinF,SinH;
  227.  
  228.     CosF=cos(F);
  229.     CosH=cos(H);
  230.     SinF=sin(F);
  231.     SinH=sin(H);
  232.     if(Z1Z2NonZeroFlag)
  233.     {
  234.         Y[1][0]=CosF*CosH-SinF*Z2+SinF*CosH*Z1;
  235.         Y[2][0]=SinF*CosH+CosF*Z2-CosF*SinH*Z1;
  236.     }
  237.     else
  238.     {
  239.         Y[1][0]=CosF*CosH;
  240.         Y[2][0]=SinF*CosH;
  241.     }
  242.     Y[3][0]=SinH;
  243. }
  244.  
  245. void Convert::SubrB(double F,double H,double Z1,double Z2)
  246. {
  247.     double CosF,CosH,SinF,SinH;
  248.  
  249.     CosF=cos(F);
  250.     CosH=cos(H);
  251.     SinF=sin(F);
  252.     SinH=sin(H);
  253.     if(Z1Z2NonZeroFlag)
  254.     {
  255.         Y[1][1]=CosF*CosH+SinF*Z2-SinF*CosH*Z1;
  256.         Y[2][1]=SinF*CosH-CosF*Z2+CosF*SinH*Z1;
  257.     }
  258.     else
  259.     {
  260.         Y[1][1]=CosF*CosH;
  261.         Y[2][1]=SinF*CosH;
  262.     }
  263.     Y[3][1]=SinH;
  264. }
  265.  
  266. void Convert::DeterminateSubr(double& W)
  267. {
  268.     W=V[1][1]*V[2][2]*V[3][3]+V[1][2]*V[2][3]*V[3][1]+V[1][3]*V[3][2]*V[2][1]-
  269.     V[1][3]*V[2][2]*V[3][1]-V[1][1]*V[3][2]*V[2][3]-V[1][2]*V[2][1]*V[3][3];
  270. }
  271.  
  272. void Convert::AngleSubr(double& F,double& H)
  273. {
  274.     double C;
  275.  
  276.     C=sqrt(Y[1][1]*Y[1][1]+Y[2][1]*Y[2][1]);
  277.     if(C==0 && Y[3][1]>0)
  278.         H=OneRev/4;
  279.     if(C==0 && Y[3][1]<0)
  280.         H=-OneRev/4;
  281.     if(C!=0)
  282.         H=atan(Y[3][1]/C);
  283.     if(C==0)
  284.     {
  285.         // F should be indeterminate: program listing is F=1000/RadToDeg;
  286.         F=0;
  287.     }
  288.     if(C!=0 && Y[1][1]==0 && Y[2][1]>0)
  289.         F=OneRev/4;
  290.     if(C!=0 && Y[1][1]==0 && Y[2][1]<0)
  291.         F=OneRev*(3/4);
  292.     if(Y[1][1]>0)
  293.         F=atan(Y[2][1]/Y[1][1]);
  294.     if(Y[1][1]<0)
  295.         F=atan(Y[2][1]/Y[1][1])+OneRev/2;
  296.     while(F>OneRev)
  297.         F-=OneRev;
  298.     while(F<0)
  299.         F+=OneRev;
  300. }
  301.  
  302. void Convert::CalcPostInitVars(void)
  303. {
  304.     Position Temp=Current;
  305.  
  306.     // find Lat by setting for 90 deg declination(NCP)
  307.     Current.Dec=OneRev/4;
  308.     Current.Ra=0;
  309.     GetAltaz();
  310.     Lat=Current.Alt;
  311.     SinLatDividedByCosLat=sin(Lat)/cos(Lat);
  312.  
  313.     // set azimuth offset
  314.     AzOff=Current.Az;
  315.  
  316.     // find longitude
  317.     Current.Alt=OneRev/4;
  318.     Current.Az=0;
  319.     GetEquat();
  320.     // ATimes::Longitude+Current.SidT.Time=Greenwich Sidereal Time;
  321.     // difference between GST and Current.Ra (==zenith) will be scope longitude
  322.     Longitude=ATimes::Longitude+Current.SidT.Time-Current.Ra;
  323.     while(Longitude<0)
  324.         Longitude+=OneRev;
  325.     while(Longitude>OneRev)
  326.         Longitude-=OneRev;
  327.  
  328.     // find hour angle offset=LST(Current.SidT.Time)-scope's meridian,
  329.     // HA=LST-HAOff-Ra, or, HA=LST-HAOff-Ra, by setting for zenith;
  330.     // +offset=scope tilted to West, -offset=scope tilted to East;
  331.     // HAOff varies from -offset to +offset(should be a small amount)
  332.     HAOff=Current.SidT.Time-Current.Ra;
  333.     while(HAOff > M_PI)
  334.         HAOff-=OneRev;
  335.     while(HAOff < -M_PI)
  336.         HAOff+=OneRev;
  337.  
  338.     // restore current coordinates
  339.     Current=Temp;
  340. }
  341.  
  342. void Convert::GetAltaz(void)
  343. {
  344.     int I,J;
  345.     double D,B,F,H;
  346.     double CosD;
  347.  
  348.     D=Current.Dec;
  349.     CosD=cos(D);
  350.     // B is CCW so HA formula backwards
  351.     B=Current.Ra-Current.SidT.Time;
  352.     X[1][1]=CosD*cos(B);
  353.     X[2][1]=CosD*sin(B);
  354.     X[3][1]=sin(D);
  355.     Y[1][1]=0;
  356.     Y[2][1]=0;
  357.     Y[3][1]=0;
  358.     for(I=1; I<=3; I++)
  359.         for(J=1; J<=3; J++)
  360.             Y[I][1]+=(R[I][J]*X[J][1]);
  361.     AngleSubr(F,H);
  362.     SubrB(F,H,Z1,Z2);
  363.     AngleSubr(F,H);
  364.     H-=Z3;
  365.     Current.Alt=H;
  366.     // convert back to CW
  367.     Current.Az=OneRev-F;
  368. }
  369.  
  370. void Convert::GetEquat(void)
  371. {
  372.     int I,J;
  373.     double F,H;
  374.  
  375.     // convert from CW to CCW az
  376.     F=OneRev-Current.Az;
  377.     H=Current.Alt+Z3;
  378.     SubrA(F,H,Z1,Z2);
  379.     X[1][1]=Y[1][0];
  380.     X[2][1]=Y[2][0];
  381.     X[3][1]=Y[3][0];
  382.     Y[1][1]=0;
  383.     Y[2][1]=0;
  384.     Y[3][1]=0;
  385.     for(I=1; I<=3; I++)
  386.         for(J=1; J<=3; J++)
  387.             Y[I][1]+=(Q[I][J]*X[J][1]);
  388.     AngleSubr(F,H);
  389.     F+=Current.SidT.Time;
  390.     Current.Ra=F;
  391.     ValidRa(Current);
  392.     Current.Dec=H;
  393. }
  394.  
  395. void Convert::CalcFieldR(void)
  396. {
  397.     double A;
  398.     double SinHA;
  399.  
  400.     HA=Current.SidT.Time-HAOff-Current.Ra;
  401.     SinHA=sin(HA);
  402.     A=SinLatDividedByCosLat*cos(Current.Dec)-sin(Current.Dec)*cos(HA);
  403.     if(A<0)
  404.         FieldR=atan(SinHA/A)+M_PI/2;
  405.     else
  406.         if(A==0)
  407.             if(SinHA<0)
  408.                 FieldR=-M_PI/2;
  409.             else
  410.                 if(SinHA==0)
  411.                     FieldR=0;
  412.                 else
  413.                     FieldR=M_PI/2;
  414.         else
  415.             FieldR=atan(SinHA/A);
  416.     if(FieldR<0)
  417.         FieldR+=OneRev;
  418. }
  419.  
  420. void Convert::DisplayFieldR(void)
  421. {
  422.     cout.setf(ios::showpoint);
  423.     cout.setf(ios::fixed,ios::floatfield);
  424.     cout<<setprecision(3)<<setw(7)<<FieldR*RadToDeg;
  425. }
  426.  
  427. void Convert::Test(void)
  428. {
  429.     clrscr();
  430.     cout<<endl<<endl<<endl<<"Test of Convert class functions:"<<endl;
  431.     InitFlag=0;
  432.     SetMountErrorsDeg(-.04,.4,-1.63);
  433.     SetCoordDeg(Current,79.172,45.998,39.9,360-39.9,39.2*SidRate/4);
  434.     Convert::InitMatrix(InitOne);
  435.     SetCoordDeg(Current,37.96,89.264,36.2,360-94.6,40.3*SidRate/4);
  436.     Convert::InitMatrix(InitTwo);
  437.     SetCoordDeg(Current,326.05,9.88,0,0,47*SidRate/4);
  438.     Convert::GetAltaz();
  439.     ShowCoordDeg(Current);
  440.     cout<<endl<<"Should be Alt: 42.16   Az: "<<(360-202.54);
  441.     SetCoordDeg(Current,71.53,17.07,0,0,62*SidRate/4);
  442.     Convert::GetAltaz();
  443.     ShowCoordDeg(Current);
  444.     cout<<endl<<"Should be Alt: 40.31   Az: "<<(360-359.98);
  445.     SetCoordDeg(Current,0,0,35.5,360-24.1,71.9*SidRate/4);
  446.     Convert::GetEquat();
  447.     ShowCoordDeg(Current);
  448.     cout<<endl<<"Should be Ra: 87.99   Dec: 32.51";
  449.     cout<<endl<<endl<<ContMsg;
  450.     getch();
  451. }
  452.  
  453.