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