help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Help-gsl] nan error in matrix-vector multiplication


From: Tim Rambo
Subject: [Help-gsl] nan error in matrix-vector multiplication
Date: Tue, 06 Jan 2009 16:58:08 -0500

The problem is this: given the two angles from a unit vector in spherical 
coordinates, I create a cartesian vector which I then multiply by the identity 
matrix. The first two elements of the vector are correct but the final element 
of the vector is nan. Any help is appreciated.
 
Thanks
 
 
 
g_rotate(double *az, double *el){
      FILE *file;
     
      double d[9];
      /*creating cartesian vector*/
      double v[]= 
{sin((*el)*M_PI/180.0)*cos((*az)*M_PI/180.0),sin((*el)*M_PI/180.0)*sin((*az)*M_PI/180.0),
 cos((*el)*M_PI/180.0)};  
      printf("%lf    %lf     %lf\n", v[0], v[1], v[2]);   
        /*reading matrix from file*/
      file = fopen("matrix.dat", "r");
      if(file==NULL){
                     printf("unable to read file\n");
                     return -1;
                     }               
  
      fscanf(file, "%lf%lf%lf%lf%lf%lf%lf%lf%lf", &d[0], &d[1], &d[2], &d[3], 
&d[4], &d[5], &d[6], &d[7], &d[8]);
  fclose(file);     
  
       printf("%lf    %lf     %lf\n %lf    %lf     %lf\n %lf    %lf     %lf ", 
d[0], d[1], d[2], d[3], d[4], d[5], d[6], d[7], d[8]);
       /*creating linear system*/      
      gsl_matrix_view R = gsl_matrix_view_array(d, 3,3);
 
      gsl_vector_view V = gsl_vector_view_array(v,3);
  
 
      gsl_vector *CV= gsl_vector_alloc(3);
     /*multiplying matrix and vector*/
      gsl_blas_dgemv(CblasNoTrans, 1.0, &R.matrix, &V.vector, 1.0, CV);
      printf("corrected vector in rectangular coordinates\n");
      gsl_vector_fprintf(stdout, CV, "%g");
       /*decomposing vecor into angles and storing the information*/ 
       double xsq=gsl_vector_get(CV, 0)*gsl_vector_get(CV, 0);
       double ysq=gsl_vector_get(CV, 1)*gsl_vector_get(CV, 1);
       double zsq=gsl_vector_get(CV, 2)*gsl_vector_get(CV, 2);       
       *az=(float)(acos( gsl_vector_get(CV, 2)/sqrt(xsq+ysq+zsq))*180.00/M_PI); 
   
       *el=(float)(atan(gsl_vector_get(CV, 1)/gsl_vector_get(CV, 
0))*180.00/M_PI);
          
        /*garbage collection*/
       gsl_vector_free(CV);
          
      }


reply via email to

[Prev in Thread] Current Thread [Next in Thread]