Calculo en coma flotante

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 :slight_smile:

El código es cortito.

/************************************************************************
 *
 ************************************************************************/
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
}

Muchas gracias.

¿Por que no vas sacando los datos por el puerto seire a medida que los calculas en el Arduino y los comparas con los calculados a mano?

La idea de @Chiva esta genial porque despues de leer coordenadas galácticas y coordenadas horizontales me he quedado a cuadros.

Los resultados de los calculos ya los estoy leyendo por el puerto serie y comparándolos con los que me da Stellarium, por ejemplo, veo una diferencia que crece a medida que me aparto de la direccion SUR.
Por eso creo que hay algún problema de redondeo en el cálculo o algo que se me pasó por alto en las fórmulas.