Conversion entre sistemas de coordenadas

[Angel Fraguela]

Saludos a todos.

Estoy haciendo un controlador para motorizar un telescopio con montura "altazimutal" para conectarlo a internet y usarlo desde programas como Stellarium o KStars.

Mi problema es que en lo más importante del asunto, que es la trigonometría esférica, debo de haber cometido algún error en las funciones de conversión entre coordenadas galácticas y coordenadas horizontales. Los resultados que obtengo no son correctos; se desvian en muchos casos varios grados.

Debe de ser algún problema con los decimales en el calculo en coma flotante. No creo que sean las formulas en si mismas porque los resultados son aproximados.

Si alguien pudiera echarle un vistazo al código de las funciones que me dan dolor de cabeza se lo agradecería eternamente.

Me da rabia haberme gastado pasta en motores y en encoders con bastante precisión para que luego los cálculos sean sólo aproximados 

El código es cortito.

Código [Seleccionar] Expandir


/************************************************************************
*
************************************************************************/
static float RAD=0.0174532925199433;  // PI/180 pero mas preciso  (no creo que se usen todos los digitos)

/*
*  J2000 a horizontales
*/
hCoord ARDEC2ALTAZ(ehCoord ardec, Fecha fecha, gCoord gc){

  float LST,RA,Dec,Lat;
  char buffer[20];
  hCoord hc;
  float AZ,HA,ALT;

  LST=LocalSiderialTime(fecha,gc);
  RA=(ardec.angulohorario.horas+(ardec.angulohorario.minutos/60.0)+(ardec.angulohorario.segundos/3600.0))*15.0;        // En grados
  Dec=(ardec.declinacion.grados+(ardec.declinacion.minutos/60.0)+(ardec.declinacion.segundos/3600.0))*RAD;             // En radians
  Lat=(gc.latitud.grados+(gc.latitud.minutos/60.0)+(gc.latitud.segundos/3600.0))*RAD;                                  // En radians

  HA=LST-RA;
  HA=(HA<0)?HA+360.0:HA;

  ALT=asin(sin(Dec)*sin(Lat)  + cos(Dec)*cos(Lat)*cos(HA*RAD));
  AZ=acos((sin(Dec)-sin(ALT)*sin(Lat))/(cos(ALT)*cos(Lat)));
  ALT=ALT/RAD;
  AZ=AZ/RAD;

  if(sin(HA*RAD) > 0)
    AZ=360.0-AZ;

  hc.altura.grados=floor(ALT);
  hc.altura.minutos=floor((ALT - floor(ALT)) * 60.0);
  hc.altura.segundos=((ALT - floor(ALT)) * 60.0 - hc.altura.minutos) * 60.0;

  hc.azimut.grados=floor(AZ);
  hc.azimut.minutos=floor((AZ - floor(AZ)) * 60.0);
  hc.azimut.segundos=((AZ - floor(AZ)) * 60.0 - hc.azimut.minutos) * 60.0;

  return hc;
}
/*
*  Horizontales a J2000
*/
ehCoord ALTAZ2ARDEC(hCoord altaz, Fecha fecha, gCoord gc){

  float LST,LAT,LON;
  ehCoord eh;
  float Dec,Ha,RA;
  int JT;
  float ALT,AZ;
 
  ALT=(altaz.altura.grados+(altaz.altura.minutos/60.0)+(altaz.altura.segundos/3600.0))*RAD;       // En radians
  AZ=(altaz.azimut.grados+(altaz.azimut.minutos/60.0)+(altaz.azimut.segundos/3600.0))*RAD;        // En radians

  LST=LocalSiderialTime(fecha,gc); // En grados

  LAT=(gc.latitud.grados+(gc.latitud.minutos/60.0)+(gc.latitud.segundos/3600.0))*RAD;                 // En radians
  LON=(gc.longitud.grados+(gc.longitud.minutos/60.0)+(gc.longitud.segundos/3600.0))*RAD;              // En radians

  Dec = asin((cos(AZ)*cos(ALT)*cos(LAT)) + (sin(ALT)*sin(LAT)));
  Ha = acos(  (sin(ALT) - (sin(Dec)*sin(LAT))) / (cos(Dec)*cos(LAT))   );
  Dec=Dec/RAD; // En grados
  Ha=Ha/RAD;   // En grados
 
  RA = LST - Ha;
  RA=RA/15.0; // En horas

  eh.declinacion.grados=floor(Dec);
  eh.declinacion.minutos=floor((Dec-floor(Dec))*60.0);
  eh.declinacion.segundos=((Dec - floor(Dec)) * 60.0 - eh.declinacion.minutos) * 60.0;

  eh.angulohorario.horas=floor(RA);
  eh.angulohorario.minutos=floor((RA-floor(RA))*60.0);
  eh.angulohorario.segundos=((RA - floor(RA)) * 60.0 - eh.angulohorario.minutos) * 60.0;

  return eh;
}
/*
*  Calculo dia juliano J2000 de una fecha determinada
*/
dJuliano diaJuliano(Fecha fc) {
  float A,Y,M;
  dJuliano dj;

  A=(14-fc.mes)/12.0;
  Y=fc.ano+4800-A;
  M=fc.mes+(12*A)-3;

  dj.dia=fc.dia+floor(((153*M)+2)/5)+floor(365*Y)+floor(Y/4)-floor(Y/100)+floor(Y/400)-32045;
  dj.fraccion=((fc.horas-12)/24.0)  +   fc.minutos/1440.0    +    fc.segundos/86400.0 ;
  dj.dia=dj.dia+floor(dj.fraccion);
  dj.fraccion=dj.fraccion-floor(dj.fraccion);
 

  return dj;
}
/*
* 
*/
float LocalSiderialTime(Fecha fc, gCoord gc){
  dJuliano dj;
  float J, lst, tiempo, lon;

  dj=diaJuliano(fc);
  tiempo=fc.horas+(fc.minutos/60.0)+(fc.segundos/3600.0);                                               // En horas. OJO: TU.
  lon=gc.longitud.grados+(gc.longitud.minutos/60.0)+(gc.longitud.segundos/3600.0);                      // En grados. Oeste negativo

  J=dj.dia-2451545.0;
  J=J+dj.fraccion;

  lst=100.46 + 0.985647 * J + lon + 15.0*tiempo ;
  lst=lst+(0.985647 * dj.fraccion);

  lst=lst-(360.0*floor(lst/360.0));

  return lst;  // En grados
}


Perdonad por el tocho.
Muchísimas gracias.

[Alejandro Quilez]

Hola y bienvenido al foro.

Como veras por el número de respuestas este es un foro de astronomia no de programacion.

De todas maneras puedes hacerle un trace para ir viendo los resultados y verificar si efectivamente son los decimales u otro el problema.


Hace años que deje el C (y la programación) y leyendo tu código apenas si recuerdo las sentencias.

Saludos.