Pages: [1]   Go Down
Author Topic: Calculo en coma flotante  (Read 810 times)
0 Members and 1 Guest are viewing this topic.
Offline Offline
Newbie
*
Karma: 0
Posts: 2
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

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  smiley

El código es cortito.
Code:
/************************************************************************
 *
 ************************************************************************/
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.
Logged

Offline Offline
Edison Member
*
Karma: 23
Posts: 1375
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

¿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?
Logged

Mercadillo electrónico. Kit iniciación a Arduino, shield LCD a color y más cosas!

0
Offline Offline
Edison Member
*
Karma: 8
Posts: 1040
Arduino rocks
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

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

Trabajando en ...

    * Control Domotico (En montaje ...)
    http://casitadomotica.blogspot.com/
 

[url=https://bitbucket.org/fmalpartida

Offline Offline
Newbie
*
Karma: 0
Posts: 2
View Profile
 Bigger Bigger  Smaller Smaller  Reset Reset

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.
Logged

Pages: [1]   Go Up
Jump to: