#include "shape.h"
 
int ecl2ast(double *rotstate, double *pos_earth, double *pos_sun, double epoch, 
	double *new_pos_earth, double *new_pos_sun) {
	
	int i ;
	double cbeta, clambda ; // Converted pole directions to radians
	double revolutions ;
	double dt ;
	double lenght_pos_earth, lenght_pos_sun ;
	double rotation_state ;
	double pos_tmp[3] ;
	double normal_pos_earth[3], normal_pos_sun[3] ;
	
	clambda = rotstate[0] ;
	cbeta = rotstate[1] ;
	revolutions = 1./rotstate[2] ;
	dt = (epoch - rotstate[3]) ;
	//printf(" --- %6.2lf  %6.2lf --- ",ra*180./PI,dec*180./PI) ;

	lenght_pos_earth = sqrt(fvector_dot(pos_earth,pos_earth)) ;
	lenght_pos_sun = sqrt(fvector_dot(pos_sun,pos_sun)) ;
	
	//printf(" %8.2lf  %8.2lf \n",pos_sun[0],lenght_pos_sun) ;
	// Normalize the vectors 
	for(i=0;i<3;i++) {
		normal_pos_earth[i] = pos_earth[i]/lenght_pos_earth ;
		normal_pos_sun[i] = pos_sun[i]/lenght_pos_sun ;
	}

	// Compute the coords of the observer and Sun in the asteroid's frame
	rotation_state = revolutions*dt*TWOPI ;
	rotation_state = modangle_radians(rotation_state) ;
	//printf("   %8.4lf  ",rotation_state) ;
	transform(clambda, cbeta, rotation_state, normal_pos_sun, new_pos_sun) ;
        //printf(" %6.2lf %6.2lf %6.2lf - ",pos_tmp[0],pos_tmp[1],pos_tmp[2]) ;
	//angpos_sun[1] = acos(pos_tmp[2]) ;
	//angpos_sun[0] = atan2(pos_tmp[1],pos_tmp[0]) ;
	transform(clambda, cbeta, rotation_state, normal_pos_earth, new_pos_earth) ;
        //printf(" %6.2lf %6.2lf %6.2lf - ",pos_tmp[0],pos_tmp[1],pos_tmp[2]) ;
	//angpos_earth[1] = acos(pos_tmp[2]) ;
	//angpos_earth[0] = atan2(pos_tmp[1],pos_tmp[0]) ;	

	return ;
}

double modangle_radians(double angle) {
  // Makes dure that the angle is between 0 and 2.*PI
  while(angle>=TWOPI) {
    angle = angle - TWOPI ;
  }
  while(angle<0.) {
    angle = angle + TWOPI ;	
  }
  return(angle) ;
}

int transform(double alpha, double beta, double gamma, double *vec0, double *vec1) {
	
	int i, j ;
	double cosa, sina, cosb, sinb, cosg, sing ;
	double tmatrix[3][3] ;

	cosa = cos(alpha) ;
	sina = sin(alpha) ;
	cosb = cos(beta) ;
	sinb = sin(beta) ;
	cosg = cos(gamma) ;
	sing = sin(gamma) ;

	tmatrix[0][0] = cosg*cosb*cosa - sing*sina ;
	tmatrix[0][1] = cosg*cosb*sina + sing*cosa ;
	tmatrix[0][2] = -cosg*sinb ;
	tmatrix[1][0] = -sing*cosb*cosa - cosg*sina ;
	tmatrix[1][1] = -sing*cosb*sina + cosg*cosa ;
	tmatrix[1][2] = sing*sinb ;
	tmatrix[2][0] = sinb*cosa ;
	tmatrix[2][1] = sinb*sina ;
	tmatrix[2][2] = cosb ;
	
	for(i=0;i<3;i++) {
		vec1[i] = 0. ;
		for(j=0;j<3;j++) vec1[i] = vec1[i] + tmatrix[i][j]*vec0[j] ;
	}
}  


