#include <stdexcept>
#include "fortran.hpp"
#include "numvec.hpp"

namespace mort
{
    using std::string;

    using std::runtime_error;

    numvec makevec( double x, double y ) 
    {
        double r[] = {x, y};
        return numvec( r, 2 );
    }

    numvec makevec( double x, double y, double z ) 
    {
        double r[] = {x, y, z};
        return numvec( r, 3 );
    }

    numvec makevec( double x, double y, double z, double s ) 
    {
        double r[] = { x, y, z, s};
        return numvec( r, 4 );
    }

    numvec makevec( double x, double y, double z, double s, double u )
    {
        double r[] = { x, y, z, s, u };
        return numvec( r, 5 );
    }

    numvec makevec( double x, double y, double z, double s, double u, double v ) 
    {
        double r[] = { x, y, z, s, u, v};
        return numvec( r, 6 );
    }
    
    double norm( const numvec& v )
    {
        if( v.size()==0 )
            throw runtime_error( "Error: in norm vector size==0" );

        return sqrt( dotpd(v,v) );
    }

    double dotpd( const numvec& a, const numvec& b )
    {
        if( a.size()!=b.size() )
            throw runtime_error( "Error: dot product two vectors difference size" );

#       ifdef MORT_DONT_USE_VALARRAY
        double accu(0);
        for (size_t n(0); n < a.size(); ++n)
            accu += a[n]*b[n];
        return accu;
#       else
        return (a*b).sum();
#       endif // MORT_DONT_USE_VALARRAY
    }    
 
    numvec cross( const numvec& v1, const numvec& v2 )    
    {
        return makevec( v1[1] * v2[2] - v2[1] * v1[2],
                        v1[2] * v2[0] - v2[2] * v1[0],
                        v1[0] * v2[1] - v2[0] * v1[1] );
    }

    numvec& normalize( numvec& vec )
    {
        vec /= norm(vec);
        return vec;
    }

    numvec normalcpy( const numvec& vec )
    {
        numvec copy( vec );
        normalize( copy );
        return copy;
    }

    
   
} // namespace mort



/*
  numvec atov( const string& str )
  {
  int curly = str.find( '{' );

  int start = (curly == (int)string::npos) ? 0 : curly + 1;
        
  std::istringstream is( str.c_str() + start );
        
  numvec v;
  double t;
  while( is >> t )
  {
  v.push_back(t);
  }
	
  return v;
  }


  numvec atov3d( const string& str )
  {
  return atov< numvec >( str );
  }
    
  vector4d atov4d( const string& str )
  {
  return atov< vector4d >( str );
  }
*/