gcc Versionsunterschied -> Problem

U

upperlimit

Jungspund
Hallo,

ich habe folgendes Problem:

Unter SuSE 10.1 (gcc 4.1.0, 64Bit, Kernel 2.6.16.13-4-default) laesst sich die Datei

Code:
//------------------------------------------------------------------------------
//
// Exercise_4_2.cpp
// 
// Purpose: 
//
//   Satellite Orbits - Models, Methods, and Applications
//   Exercise 4-2: Gauss-Jackson 4th-order predictor 
//
// Notes:
//
//   This software is protected by national and international copyright. 
//   Any unauthorized use, reproduction or modificaton is unlawful and 
//   will be prosecuted. Commercial and non-private application of the 
//   software in any form is strictly prohibited unless otherwise granted
//   by the authors.
//
//   The code is provided without any warranty; without even the implied 
//   warranty of merchantibility or fitness for a particular purpose.
//
// Last modified:
//
//   2000/03/04  OMO  Final version (1st edition)
//   2001/08/17  OMO  Minor upgrades
//
// (c) 1999-2001  O. Montenbruck, E. Gill
//
//------------------------------------------------------------------------------

#include <cmath>
#include <iostream>
#include <iomanip>

#include "SAT_Const.h"
#include "SAT_Kepler.h"
#include "SAT_VecMat.h"

using namespace std;


//------------------------------------------------------------------------------
//
// GJ4P class 
//
//------------------------------------------------------------------------------

// Function prototype for second order differential equations
// void f (double t, const Vector& r, const Vector& v, Vector& a)

typedef void (*GJ4Pfunct)(
  double        t,     // Independent variable
  const Vector& r,     // Position vector 
  const Vector& v,     // Velocity vector  r'=v
  Vector&       a,     // Acceleration     r''=a=f(t,r,v)
  void*         pAux   // Pointer to auxiliary data used within f
);


// Specification

class GJ4P
{
  public:

    // Constructor
    GJ4P (
      GJ4Pfunct   f_,        // Differential equation
      int         n_eqn_,    // Dimension
      void*       pAux_      // Pointer to auxiliary data
      )
    : n_eqn(n_eqn_), f(f_), pAux(pAux_)
    {};
    
    // Initialization step
    void Init (
      double        t_0,     // Initial value of independent variable
      const Vector& r_0,     // Initial value r_0=r(t_0)
      const Vector& v_0,     // Initial value v_0=dr/dt(t_0)
      double        h_       // Step size
    );

    // Integration step
    void Step (         
      double&  t,            // Independent variable; updated by t+h
      Vector&  r,            // Value of r(t); updated by r(t+h)
      Vector&  v             // Value of v(t)=dr/dt(t); updated by v(t+h)
    );

  private:

    // 4th order Runge-Kutta step
    void RK4 (         
      double&  t,            // Independent variable; updated by t+h
      Vector&  r,            // Value of r(t); updated by r(t+h)
      Vector&  v,            // Value of v(t)=dr/dt(t); updated by v(t+h)
      double   h             // Step size
    );

    // Elements
    int         n_eqn;       // Dimension
    GJ4Pfunct   f;           // Differential equation
    double      h;           // Step size
    void*       pAux;        // Pointer to auxiliary data requird by f
    Vector      S2,S1;       // First and second sum of acceleration
    Vector      D[4];        // Backward differences of acceleration at t
    Vector      d[4];        // Backward differences of acceleration at t+h
    Vector      r_p,v_p;     // Predictor

};


//
// 4th order Runge-Kutta step for 2nd order differential equation
//

void GJ4P::RK4 (double&  t, Vector&  r, Vector&  v, double h )
{
  Vector v_1, v_2, v_3, v_4;
  Vector a_1, a_2, a_3, a_4;

  v_1 = v;              f( t      , r            , v_1, a_1, pAux ); 
  v_2 = v+(h/2.0)*a_1;  f( t+h/2.0, r+(h/2.0)*v_1, v_2, a_2, pAux );
  v_3 = v+(h/2.0)*a_2;  f( t+h/2.0, r+(h/2.0)*v_2, v_3, a_3, pAux ); 
  v_4 = v+h*a_3;        f( t+h    , r+h*v_3      , v_4, a_4, pAux );
  
  t = t + h;
  r = r + (h/6.0)*( v_1 + 2.0*v_2 + 2.0*v_3 + v_4 );
  v = v + (h/6.0)*( a_1 + 2.0*a_2 + 2.0*a_3 + a_4 );

};


//
// Initialization of backwards differences from initial conditions
//

void GJ4P::Init(double t_0, const Vector& r_0, const Vector& v_0, double h_)
{
  // Order of method
  
  const int m = 4;

  // Coefficients gamma/delta of 1st/2nd order Moulton/Cowell corrector method

  const double gc[m+1] = {+1.0, -1/2.0, -1/12.0, -1/24.0, -19/720.0 };
  const double dc[m+2] = {+1.0,   -1.0, +1/12.0,     0.0,  -1/240.0, -1/240.0 };

  int       i,j;
  double    t = t_0;
  Vector    r = r_0;
  Vector    v = v_0;

  // Save step size  

  h = h_;     

  // Create table of accelerations at past times t-3h, t-2h, and t-h using
  // RK4 steps

  f(t,r,v,D[0],pAux);     // D[i]=a(t-ih)
  for (i=1;i<=m-1;i++) {
    RK4(t,r,v,-h);  f(t,r,v,D[i],pAux);   
  };

  // Compute backwards differences
  
  for (i=1;i<=m-1;i++) 
    for (j=m-1;j>=i;j--) D[j] = D[j-1]-D[j];

  // Initialize backwards sums using 4th order GJ corrector

  S1 = v_0/h;              for (i=1;i<=m  ;i++) S1 -= gc[i]*D[i-1];
  S2 = r_0/(h*h)-dc[1]*S1; for (i=2;i<=m+1;i++) S2 -= dc[i]*D[i-2];

};


//
// Step from t to t+h
//

void GJ4P::Step (double& t, Vector& r, Vector& v) 
{
  // Order of method
  
  const int m = 4;  

  // Coefficients gamma/delta of 1st/2nd order Bashforth/Stoermr predictor
  
  const double gp[m+1] = {+1.0, +1/2.0, +5/12.0,  +3/8.0, +251/720.0 };
  const double dp[m+2] = {+1.0,    0.0, +1/12.0, +1/12.0,  +19/240.0,  +3/40.0 };

  int i;
  
  // 4th order predictor

  r_p = dp[0]*S2; for(i=2;i<=m+1;i++) r_p += dp[i]*D[i-2]; r_p = (h*h)*r_p;
  v_p = gp[0]*S1; for(i=1;i<=m  ;i++) v_p += gp[i]*D[i-1]; v_p =     h*v_p;

  // Update backwards difference table

  f ( t+h, r_p,v_p, d[0], pAux );               // Acceleration at t+h
  for (i=1;i<=m-1;i++) d[i]=d[i-1]-D[i-1];      // New differences at t+h
  for (i=0;i<=m-1;i++) D[i]=d[i];               // Update differences 
  S1 += d[0];  S2 += S1;                        // Update sums

  // Update independent variable and solution

  t = t + h;
  r = r_p;
  v = v_p;

};


//------------------------------------------------------------------------------
//
// f_Kep3D
//
// Purpose:
// 
//   Computes the second time derivative of the position vector for the 
//   normalized (GM=1) Kepler's problem in three dimensions
//
// Note:
//
//   pAux is expected to point to an integer variable that will be incremented
//   by one on each call of f_Kep3D
//
//------------------------------------------------------------------------------

void f_Kep3D ( double t, const Vector& r, const Vector& v, 
               Vector& a, void* pAux )
{
  
  // Pointer to auxiliary integer variable used as function call counter
  int* pCalls = static_cast<int*>(pAux);
  
  // 2nd order derivative d^2(r)/dt^2 of the position vector
  a = -r/(pow(Norm(r),3));
  
  // Increment function call count
  (*pCalls)++;

};


//------------------------------------------------------------------------------
//
// Main program
//
//------------------------------------------------------------------------------

int main() {

  // Constants
  
  const double  GM    = 1.0;                   // Gravitational coefficient
  const double  e     = 0.1;                   // Eccentricity
  const double  t_end = 20.0;                  // End time
  const Vector  Kep(1.0,e,0.0,0.0,0.0,0.0);    // (a,e,i,Omega,omega,M)
  const Vector  y_ref = State(GM,Kep,t_end);   // Reference solution

  const int     Steps[] = { 100, 300, 600, 1000, 1500, 2000, 3000, 4000 };
  
  // Variables
  
  int     nCalls;                              // Function call count
  int     iCase;
  double  t,h;                                 // Time and step size
  Vector  y(6);                                // State vector
  Vector  r(3);
  Vector  v(3);

  GJ4P    Orbit(f_Kep3D,3,&nCalls);            // Object for integrating the
                                               // 2nd order diff. equation
                                               // defined by f_Kep3D using the
                                               // 4th-order GJ predictor
  
  // Header 

  cout << "Exercise 4-2: Gauss-Jackson 4th-order predictor" 
       << endl << endl;
  cout << "  Problem D1 (e=0.1)" << endl << endl;
  cout << "  N_fnc   Accuracy   Digits " << endl;
    
  // Loop over test cases

  for (iCase=0; iCase<8; iCase++) {
  
    // Step size
    h = t_end/Steps[iCase];

    // Initial values
    t = 0.0;
    r = Vector( 1.0-e, 0.0, 0.0 );
    v = Vector( 0.0, sqrt((1+e)/(1-e)), 0.0 );
    nCalls = 0;

    // Integration from t=t to t=t_end
    Orbit.Init ( t, r, v, h );
    for (int i=1; i<=Steps[iCase]; i++)  
      Orbit.Step( t, r, v );
    y = Stack(r,v);

    // Output
    cout << fixed  << setw(6) << nCalls
         << scientific << setprecision(3) << setw(13)
         << Norm(y-y_ref) 
         << fixed << setprecision(2) << setw(7)
         << -log10(Norm(y-y_ref)) << endl;
  
  };

  return 0;

}

( http://www.fjtakalab.ynu.ac.jp/fjtaka_lec/VLBI/Source_Astro/Source/Exercise_4_2.cpp )

mit den Header-Files

http://www.fjtakalab.ynu.ac.jp/fjtaka_lec/VLBI/Source_Astro/Source/SAT_Const.h
http://www.fjtakalab.ynu.ac.jp/fjtaka_lec/VLBI/Source_Astro/Source/SAT_Kepler.h
und
http://www.fjtakalab.ynu.ac.jp/fjtaka_lec/VLBI/Source_Astro/Source/SAT_VecMat.h

nicht kompilieren.


Der Befehl

Code:
$gcc -I. -c Exercise_4_2.cpp

bzw.
Code:
$g++ -I. -c Exercise_4_2.cpp

bringt die Fehlermeldung:

Code:
SAT_VecMat.h:106: error: ‘Matrix’ does not name a type
SAT_VecMat.h:109: error: expected ‘,’ or ‘...’ before ‘&’ token
SAT_VecMat.h:109: error: ISO C++ forbids declaration of ‘Matrix’ with no type
SAT_VecMat.h:109: error: ‘Vector operator*(int)’ must have an argument of class or enumerated type
SAT_VecMat.h:110: error: expected ‘,’ or ‘...’ before ‘&’ token
SAT_VecMat.h:110: error: ISO C++ forbids declaration of ‘Matrix’ with no type
SAT_VecMat.h:113: error: ‘Matrix’ does not name a type

Interessanterweise wird dieselbe Aufgabe mit dem gcc 3.3.3 (SuSE 9.0, 32Bit, Kernel 2.4.x-default) anstandslos erledigt.

Meine Frage lautet also:
Wie kann ich es erreichen, dass die C++-Quelle auch unter gcc 4.1.0 compiliert wird ? (Um einen Unterschied zw. 32- und 64Bit-Architektur scheint es nicht zu gehen, da dasselbe Problem auch auf einer 32Bit-Maschine auftritt, allerdings auch mit gcc 4.1.0.)

Jeder Vorschlag ist willkommen.

Gruesse,
upperlimit
 
Hallo
Wenn der Code nunmal für den gcc 3.3.3 geschrieben wurde, dann musst du ihn entweder selbst bearbeiten - der neueren Version anpassen - , oder eben mit dem erforderlichen Compiler compilieren.
Dafür gibt es ja in der Regel die Hinweise, welche Version benötigt wird. Steht meist in der README.

Gruß Wolfgang
 
Wolfgang schrieb:
Hallo
Wenn der Code nunmal für den gcc 3.3.3 geschrieben wurde, dann musst du ihn entweder selbst bearbeiten - der neueren Version anpassen - , oder eben mit dem erforderlichen Compiler compilieren.
Dafür gibt es ja in der Regel die Hinweise, welche Version benötigt wird. Steht meist in der README.

Gruß Wolfgang
Ich wuerde mal sagen, dass der Code nichts Versionsspezifisches enthaelt. gcc-4.0 ist einfach strenger und enger am Standard. In der Datei SAT_VecMat.h wird die Klasse Matrix innerhalb der Klasse Vector benutzt, ohne sie vorher zu deklarieren. Schreibe einfach ueber Zeile 44, also vor der Deklaration der Klasse Vector, ein 'class Matrix;' hin, dann sollte zumindest diese Fehlermeldung verschwinden. Es wundert mich, dass gcc-3.3 das hat durchgehen lassen.
 
Hallo,

nach einer Parallelinstallation von gcc 4.1.0 und gcc 3.3.6 (SuSE 10.0, 64Bit, Kernel 2.6.16.13-4-default) habe ich es nun hinbekommen, dass die Datei mit dem gcc 3.3.6 kompiliert werden kann.

Dein Tipp, rikola, mit dem Einfuegen von 'class Matrix;' in die SAT_VecMat.h fuehrt allerdings zu anderen Fehler in der Datei SAT_RefSys1.cpp

Code:
//------------------------------------------------------------------------------
//
// SAT_RefSys1.cpp
// 
// Purpose: 
//
//   Transformations betweeen celestial and terrestrial reference systems
//
// Notes:
//
//   This software is protected by national and international copyright. 
//   Any unauthorized use, reproduction or modificaton is unlawful and 
//   will be prosecuted. Commercial and non-private application of the 
//   software in any form is strictly prohibited unless otherwise granted
//   by the authors.
//
//   The code is provided without any warranty; without even the implied 
//   warranty of merchantibility or fitness for a particular purpose.
//
// Last modified:
//
//   2000/03/04  OMO  Final version (1st edition)
//   
// (c) 1999-2000  O. Montenbruck, E. Gill
//
//------------------------------------------------------------------------------

#include <cmath>

#ifdef __GNUC__   // GNU C++ adaptation
#include <float.h>
#else             // Standard C++ version
#include <limits>
#endif

#include "SAT_Const.h"
#include "SAT_RefSys1.h"
#include "SAT_VecMat.h"

using namespace std;



//
// Local declarations
//

namespace { 
  // Machine accuracy
  #ifdef __GNUC__   // GNU C++ adaptation
  const double eps_mach = DBL_EPSILON;
  #else             // Standard C++ version
  const double eps_mach = numeric_limits<double>::epsilon();
  #endif
  // Fractional part of a number (y=x-[x])
  double Frac (double x) { return x-floor(x); };
  // x mod y
  double Modulo (double x, double y) { return y*Frac(x/y); }
}            


//------------------------------------------------------------------------------
//
// IERS (Class implementation)
//
// Purpose:
//
//   Management of IERS time and polar motion data
//  
//   p.157-169
//   5.1 Time
//------------------------------------------------------------------------------

// Class constants

const double IERS::TT_TAI  = +32.184; // TT-TAI time difference [s] p.162, 5.3
const double IERS::GPS_TAI = -19.0;   // GPS-TAI time difference [s] p.162, 5.4

// Default values of Earth rotation parameters

double IERS::UT1_UTC_ = 0.0;          // UT1-UTC time difference [s]
double IERS::UTC_TAI_ = 0.0;          // UTC-TAI time difference [s] 
double IERS::x_pole_  = 0.0;          // Pole coordinate [rad]
double IERS::y_pole_  = 0.0;          // Pole coordinate [rad]

// Setting of IERS Earth rotation parameters
// (UT1-UTC [s],UTC-TAI [s], x ["], y ["])
// cf. Fig.5.4, p.166, Buuletin B of the IERS 

void IERS::Set(double UT1_UTC, double UTC_TAI,
               double x_pole, double y_pole) 
{
   UT1_UTC_ = UT1_UTC;
   UTC_TAI_ = UTC_TAI;
   x_pole_  = x_pole/Arcs;
   y_pole_  = y_pole/Arcs;
};

// UT1-UTC time difference [s]

double IERS::UT1_UTC(double Mjd_UTC) { return UT1_UTC_; };
     
// UTC-TAI time difference [s]

double IERS::UTC_TAI(double Mjd_UTC) { return UTC_TAI_; };

// TT-UTC time difference [s]

double IERS::TT_UTC(double Mjd_UTC) { return TT_TAI-UTC_TAI_; };

// GPS-UTC time difference [s]

double IERS::GPS_UTC(double Mjd_UTC) { return GPS_TAI-UTC_TAI_; };

// Pole coordinate [rad]

double IERS::x_pole(double Mjd_UTC) { return x_pole_; };

// Pole coordinate [rad]

double IERS::y_pole(double Mjd_UTC) { return y_pole_; };



//------------------------------------------------------------------------------
//
// MeanObliquity
//
// Purpose:
//
//   Computes the mean obliquity of the ecliptic
//
// Input/Output:
//
//   Mjd_TT    Modified Julian Date (Terrestrial Time)
//   <return>  Mean obliquity of the ecliptic
//
//  p.176
//  5.3.2 Coordinate Changes due to Precession
//  5.42
//------------------------------------------------------------------------------

double MeanObliquity (double Mjd_TT)
{
  const double T = (Mjd_TT-MJD_J2000)/36525.0;

  return 
    Rad *( 23.43929111-(46.8150+(0.00059-0.001813*T)*T)*T/3600.0 );
}


//------------------------------------------------------------------------------
//
// EclMatrix
//
// Purpose:
//
//   Transformation of equatorial to ecliptical coordinates
//
// Input/Output:
//
//   Mjd_TT    Modified Julian Date (Terrestrial Time)
//   <return>  Transformation matrix
//
//  p.170
//  5.2 Celestial and Terrestrial Reference Systems
//  5.21
//------------------------------------------------------------------------------

Matrix EclMatrix (double Mjd_TT)
{
  return R_x ( MeanObliquity(Mjd_TT) );
}


//------------------------------------------------------------------------------
//
// PrecMatrix
//
// Purpose:
//
//   Precession transformation of equatorial coordinates
//
// Input/Output:
//
//   Mjd_1     Epoch given (Modified Julian Date TT)
//   MjD_2     Epoch to precess to (Modified Julian Date TT)
//   <return>  Precession transformation matrix
//
//  p.174-177
//  5.3.2 Coordinate Changes due to Precession
//  5.44,5.46,5.47
//------------------------------------------------------------------------------

Matrix PrecMatrix (double Mjd_1, double Mjd_2)
{

  // Constants

  const double T  = (Mjd_1-MJD_J2000)/36525.0;
  const double dT = (Mjd_2-Mjd_1)/36525.0;
  
  // Variables

  double zeta,z,theta;

  // Precession angles
  
  zeta  =  ( (2306.2181+(1.39656-0.000139*T)*T)+
                ((0.30188-0.000344*T)+0.017998*dT)*dT )*dT/Arcs;
  z     =  zeta + ( (0.79280+0.000411*T)+0.000205*dT)*dT*dT/Arcs;
  theta =  ( (2004.3109-(0.85330+0.000217*T)*T)-
                ((0.42665+0.000217*T)+0.041833*dT)*dT )*dT/Arcs;

  // Precession matrix
  
  return R_z(-z) * R_y(theta) * R_z(-zeta);

}    


//------------------------------------------------------------------------------
//
// NutAngles 
//
// Purpose:
//
//   Nutation in longitude and obliquity
//
// Input/Output:
//
//   Mjd_TT    Modified Julian Date (Terrestrial Time)
//   <return>  Nutation matrix
//
//  p.178-181
//  5.3.3 Nutation
//  Table 5.2, p.179, the IAU 1980 nutation theory, 5.58,5.59,5.60
//------------------------------------------------------------------------------

void NutAngles (double Mjd_TT, double& dpsi, double& deps)
{

  // Constants

  const double T  = (Mjd_TT-MJD_J2000)/36525.0;
  const double T2 = T*T;
  const double T3 = T2*T;
  const double rev = 360.0*3600.0;  // arcsec/revolution

  const int  N_coeff = 106;
  const long C[N_coeff][9] =
  {
   //
   // l  l' F  D Om    dpsi    *T     deps     *T       #
   //
    {  0, 0, 0, 0, 1,-1719960,-1742,  920250,   89 },   //   1
    {  0, 0, 0, 0, 2,   20620,    2,   -8950,    5 },   //   2
    { -2, 0, 2, 0, 1,     460,    0,    -240,    0 },   //   3
    {  2, 0,-2, 0, 0,     110,    0,       0,    0 },   //   4
    { -2, 0, 2, 0, 2,     -30,    0,      10,    0 },   //   5
    {  1,-1, 0,-1, 0,     -30,    0,       0,    0 },   //   6
    {  0,-2, 2,-2, 1,     -20,    0,      10,    0 },   //   7
    {  2, 0,-2, 0, 1,      10,    0,       0,    0 },   //   8
    {  0, 0, 2,-2, 2, -131870,  -16,   57360,  -31 },   //   9
    {  0, 1, 0, 0, 0,   14260,  -34,     540,   -1 },   //  10
    {  0, 1, 2,-2, 2,   -5170,   12,    2240,   -6 },   //  11
    {  0,-1, 2,-2, 2,    2170,   -5,    -950,    3 },   //  12
    {  0, 0, 2,-2, 1,    1290,    1,    -700,    0 },   //  13
    {  2, 0, 0,-2, 0,     480,    0,      10,    0 },   //  14
    {  0, 0, 2,-2, 0,    -220,    0,       0,    0 },   //  15
    {  0, 2, 0, 0, 0,     170,   -1,       0,    0 },   //  16
    {  0, 1, 0, 0, 1,    -150,    0,      90,    0 },   //  17
    {  0, 2, 2,-2, 2,    -160,    1,      70,    0 },   //  18
    {  0,-1, 0, 0, 1,    -120,    0,      60,    0 },   //  19
    { -2, 0, 0, 2, 1,     -60,    0,      30,    0 },   //  20
    {  0,-1, 2,-2, 1,     -50,    0,      30,    0 },   //  21
    {  2, 0, 0,-2, 1,      40,    0,     -20,    0 },   //  22
    {  0, 1, 2,-2, 1,      40,    0,     -20,    0 },   //  23
    {  1, 0, 0,-1, 0,     -40,    0,       0,    0 },   //  24
    {  2, 1, 0,-2, 0,      10,    0,       0,    0 },   //  25
    {  0, 0,-2, 2, 1,      10,    0,       0,    0 },   //  26
    {  0, 1,-2, 2, 0,     -10,    0,       0,    0 },   //  27
    {  0, 1, 0, 0, 2,      10,    0,       0,    0 },   //  28
    { -1, 0, 0, 1, 1,      10,    0,       0,    0 },   //  29
    {  0, 1, 2,-2, 0,     -10,    0,       0,    0 },   //  30
    {  0, 0, 2, 0, 2,  -22740,   -2,    9770,   -5 },   //  31
    {  1, 0, 0, 0, 0,    7120,    1,     -70,    0 },   //  32
    {  0, 0, 2, 0, 1,   -3860,   -4,    2000,    0 },   //  33
    {  1, 0, 2, 0, 2,   -3010,    0,    1290,   -1 },   //  34
    {  1, 0, 0,-2, 0,   -1580,    0,     -10,    0 },   //  35
    { -1, 0, 2, 0, 2,    1230,    0,    -530,    0 },   //  36
    {  0, 0, 0, 2, 0,     630,    0,     -20,    0 },   //  37
    {  1, 0, 0, 0, 1,     630,    1,    -330,    0 },   //  38
    { -1, 0, 0, 0, 1,    -580,   -1,     320,    0 },   //  39
    { -1, 0, 2, 2, 2,    -590,    0,     260,    0 },   //  40
    {  1, 0, 2, 0, 1,    -510,    0,     270,    0 },   //  41
    {  0, 0, 2, 2, 2,    -380,    0,     160,    0 },   //  42
    {  2, 0, 0, 0, 0,     290,    0,     -10,    0 },   //  43
    {  1, 0, 2,-2, 2,     290,    0,    -120,    0 },   //  44
    {  2, 0, 2, 0, 2,    -310,    0,     130,    0 },   //  45
    {  0, 0, 2, 0, 0,     260,    0,     -10,    0 },   //  46
    { -1, 0, 2, 0, 1,     210,    0,    -100,    0 },   //  47
    { -1, 0, 0, 2, 1,     160,    0,     -80,    0 },   //  48
    {  1, 0, 0,-2, 1,    -130,    0,      70,    0 },   //  49
    { -1, 0, 2, 2, 1,    -100,    0,      50,    0 },   //  50
    {  1, 1, 0,-2, 0,     -70,    0,       0,    0 },   //  51
    {  0, 1, 2, 0, 2,      70,    0,     -30,    0 },   //  52
    {  0,-1, 2, 0, 2,     -70,    0,      30,    0 },   //  53
    {  1, 0, 2, 2, 2,     -80,    0,      30,    0 },   //  54
    {  1, 0, 0, 2, 0,      60,    0,       0,    0 },   //  55
    {  2, 0, 2,-2, 2,      60,    0,     -30,    0 },   //  56
    {  0, 0, 0, 2, 1,     -60,    0,      30,    0 },   //  57
    {  0, 0, 2, 2, 1,     -70,    0,      30,    0 },   //  58
    {  1, 0, 2,-2, 1,      60,    0,     -30,    0 },   //  59
    {  0, 0, 0,-2, 1,     -50,    0,      30,    0 },   //  60
    {  1,-1, 0, 0, 0,      50,    0,       0,    0 },   //  61
    {  2, 0, 2, 0, 1,     -50,    0,      30,    0 },   //  62
    {  0, 1, 0,-2, 0,     -40,    0,       0,    0 },   //  63
    {  1, 0,-2, 0, 0,      40,    0,       0,    0 },   //  64
    {  0, 0, 0, 1, 0,     -40,    0,       0,    0 },   //  65
    {  1, 1, 0, 0, 0,     -30,    0,       0,    0 },   //  66
    {  1, 0, 2, 0, 0,      30,    0,       0,    0 },   //  67
    {  1,-1, 2, 0, 2,     -30,    0,      10,    0 },   //  68
    { -1,-1, 2, 2, 2,     -30,    0,      10,    0 },   //  69
    { -2, 0, 0, 0, 1,     -20,    0,      10,    0 },   //  70
    {  3, 0, 2, 0, 2,     -30,    0,      10,    0 },   //  71
    {  0,-1, 2, 2, 2,     -30,    0,      10,    0 },   //  72
    {  1, 1, 2, 0, 2,      20,    0,     -10,    0 },   //  73
    { -1, 0, 2,-2, 1,     -20,    0,      10,    0 },   //  74
    {  2, 0, 0, 0, 1,      20,    0,     -10,    0 },   //  75
    {  1, 0, 0, 0, 2,     -20,    0,      10,    0 },   //  76
    {  3, 0, 0, 0, 0,      20,    0,       0,    0 },   //  77
    {  0, 0, 2, 1, 2,      20,    0,     -10,    0 },   //  78
    { -1, 0, 0, 0, 2,      10,    0,     -10,    0 },   //  79
    {  1, 0, 0,-4, 0,     -10,    0,       0,    0 },   //  80
    { -2, 0, 2, 2, 2,      10,    0,     -10,    0 },   //  81
    { -1, 0, 2, 4, 2,     -20,    0,      10,    0 },   //  82
    {  2, 0, 0,-4, 0,     -10,    0,       0,    0 },   //  83
    {  1, 1, 2,-2, 2,      10,    0,     -10,    0 },   //  84
    {  1, 0, 2, 2, 1,     -10,    0,      10,    0 },   //  85
    { -2, 0, 2, 4, 2,     -10,    0,      10,    0 },   //  86
    { -1, 0, 4, 0, 2,      10,    0,       0,    0 },   //  87
    {  1,-1, 0,-2, 0,      10,    0,       0,    0 },   //  88
    {  2, 0, 2,-2, 1,      10,    0,     -10,    0 },   //  89
    {  2, 0, 2, 2, 2,     -10,    0,       0,    0 },   //  90
    {  1, 0, 0, 2, 1,     -10,    0,       0,    0 },   //  91
    {  0, 0, 4,-2, 2,      10,    0,       0,    0 },   //  92
    {  3, 0, 2,-2, 2,      10,    0,       0,    0 },   //  93
    {  1, 0, 2,-2, 0,     -10,    0,       0,    0 },   //  94
    {  0, 1, 2, 0, 1,      10,    0,       0,    0 },   //  95
    { -1,-1, 0, 2, 1,      10,    0,       0,    0 },   //  96
    {  0, 0,-2, 0, 1,     -10,    0,       0,    0 },   //  97
    {  0, 0, 2,-1, 2,     -10,    0,       0,    0 },   //  98
    {  0, 1, 0, 2, 0,     -10,    0,       0,    0 },   //  99
    {  1, 0,-2,-2, 0,     -10,    0,       0,    0 },   // 100
    {  0,-1, 2, 0, 1,     -10,    0,       0,    0 },   // 101
    {  1, 1, 0,-2, 1,     -10,    0,       0,    0 },   // 102
    {  1, 0,-2, 2, 0,     -10,    0,       0,    0 },   // 103
    {  2, 0, 0, 2, 0,      10,    0,       0,    0 },   // 104
    {  0, 0, 2, 4, 2,     -10,    0,       0,    0 },   // 105
    {  0, 1, 0, 1, 0,      10,    0,       0,    0 }    // 106
   };

  // Variables

  double  l, lp, F, D, Om;
  double  arg;
  

  // Mean arguments of luni-solar motion
  //
  //   l   mean anomaly of the Moon
  //   l'  mean anomaly of the Sun
  //   F   mean argument of latitude
  //   D   mean longitude elongation of the Moon from the Sun 
  //   Om  mean longitude of the ascending node
  
  l  = Modulo (  485866.733 + (1325.0*rev +  715922.633)*T    
                                 + 31.310*T2 + 0.064*T3, rev );
  lp = Modulo ( 1287099.804 + (  99.0*rev + 1292581.224)*T             
                                 -  0.577*T2 - 0.012*T3, rev );
  F  = Modulo (  335778.877 + (1342.0*rev +  295263.137)*T    
                                 - 13.257*T2 + 0.011*T3, rev );
  D  = Modulo ( 1072261.307 + (1236.0*rev + 1105601.328)*T    
                                 -  6.891*T2 + 0.019*T3, rev );
  Om = Modulo (  450160.280 - (   5.0*rev +  482890.539)*T    
                                 +  7.455*T2 + 0.008*T3, rev );

  // Nutation in longitude and obliquity [rad]

  deps = dpsi = 0.0;
  for (int i=0; i<N_coeff; i++) {
    arg  =  ( C[i][0]*l+C[i][1]*lp+C[i][2]*F+C[i][3]*D+C[i][4]*Om ) / Arcs;
    dpsi += ( C[i][5]+C[i][6]*T ) * sin(arg);
    deps += ( C[i][7]+C[i][8]*T ) * cos(arg);
  };
      
  dpsi = 1.0E-5 * dpsi/Arcs;
  deps = 1.0E-5 * deps/Arcs;

}


//------------------------------------------------------------------------------
//
// NutMatrix 
//
// Purpose:
//
//   Transformation from mean to true equator and equinox
//
// Input/Output:
//
//   Mjd_TT    Modified Julian Date (Terrestrial Time)
//   <return>  Nutation matrix
//
//  p.178-181
//  5.3.3 Nutation
//  5.62
//------------------------------------------------------------------------------

Matrix NutMatrix (double Mjd_TT)
{

  double dpsi, deps, eps;
  
  // Mean obliquity of the ecliptic

  eps = MeanObliquity(Mjd_TT);

  // Nutation in longitude and obliquity

  NutAngles (Mjd_TT, dpsi,deps);

  // Transformation from mean to true equator and equinox

  return  R_x(-eps-deps)*R_z(-dpsi)*R_x(+eps);

}


//------------------------------------------------------------------------------
//
// NutMatrixSimple
//
// Purpose:
//
//   Transformation from mean to true equator and equinox (low precision)
//
// Input/Output:
//
//   Mjd_TT    Modified Julian Date (Terrestrial Time)
//   <return>  Nutation matrix
//
//------------------------------------------------------------------------------

Matrix NutMatrixSimple (double Mjd_TT)
{

  // Constants

  const double T  = (Mjd_TT-MJD_J2000)/36525.0;

  // Variables

  double  ls, D, F, N;
  double  eps, dpsi, deps;

  // Mean arguments of luni-solar motion
  
  ls = pi2*Frac(0.993133+  99.997306*T);   // mean anomaly Sun          
  D  = pi2*Frac(0.827362+1236.853087*T);   // diff. longitude Moon-Sun  
  F  = pi2*Frac(0.259089+1342.227826*T);   // mean argument of latitude 
  N  = pi2*Frac(0.347346-   5.372447*T);   // longit. ascending node    

  // Nutation angles
  
  dpsi = ( -17.200*sin(N)   - 1.319*sin(2*(F-D+N)) - 0.227*sin(2*(F+N))
           + 0.206*sin(2*N) + 0.143*sin(ls) ) / Arcs;
  deps = ( + 9.203*cos(N)   + 0.574*cos(2*(F-D+N)) + 0.098*cos(2*(F+N))
           - 0.090*cos(2*N)                 ) / Arcs;

  // Mean obliquity of the ecliptic

  eps  = 0.4090928-2.2696E-4*T;    

  return  R_x(-eps-deps)*R_z(-dpsi)*R_x(+eps);

}


//------------------------------------------------------------------------------
//
// EqnEquinox 
//
// Purpose:
//
//   Computation of the equation of the equinoxes
//
// Input/Output:
//
//   Mjd_TT    Modified Julian Date (Terrestrial Time)
//   <return>  Equation of the equinoxes
//
// Notes:
//
//   The equation of the equinoxes dpsi*cos(eps) is the right ascension of the 
//   mean equinox referred to the true equator and equinox and is equal to the 
//   difference between apparent and mean sidereal time.
//
//  p.181
//  5.4.1 Rotation About the Celestial Ephemeris Pole
//  5.66
//------------------------------------------------------------------------------

double EqnEquinox (double Mjd_TT)
{
  double dpsi, deps;              // Nutation angles

  // Nutation in longitude and obliquity 
      
  NutAngles (Mjd_TT, dpsi,deps );

  // Equation of the equinoxes

  return  dpsi * cos ( MeanObliquity(Mjd_TT) );

};


//------------------------------------------------------------------------------
//
// GMST
//
// Purpose:
//
//   Greenwich Mean Sidereal Time
//
// Input/Output:
//
//   Mjd_UT1   Modified Julian Date UT1
//   <return>  GMST in [rad]
//
//  p.165-169
//  Sidereal Time and Universal Time
//  5.15 -5.20
//------------------------------------------------------------------------------

double GMST (double Mjd_UT1)
{

  // Constants

  const double Secs = 86400.0;        // Seconds per day

  // Variables

  double Mjd_0,UT1,T_0,T,gmst;

  // Mean Sidereal Time
  
  Mjd_0 = floor(Mjd_UT1);
  UT1   = Secs*(Mjd_UT1-Mjd_0);          // [s]
  T_0   = (Mjd_0  -MJD_J2000)/36525.0; 
  T     = (Mjd_UT1-MJD_J2000)/36525.0; 

  gmst  = 24110.54841 + 8640184.812866*T_0 + 1.002737909350795*UT1
          + (0.093104-6.2e-6*T)*T*T; // [s]

  return  pi2*Frac(gmst/Secs);       // [rad], 0..2pi

}


//------------------------------------------------------------------------------
//
// GAST
//
// Purpose:
//
//   Greenwich Apparent Sidereal Time
//
// Input/Output:
//
//   Mjd_UT1   Modified Julian Date UT1
//   <return>  GMST in [rad]
//
//  p.181
//  5.4.1 Rotation About the Celestial Ephemeris Pole
//  5.66
//------------------------------------------------------------------------------

double GAST (double Mjd_UT1)
{
  return Modulo ( GMST(Mjd_UT1) + EqnEquinox(Mjd_UT1), pi2 );
}


//------------------------------------------------------------------------------
//
// GHAMatrix
//
// Purpose:
//
//   Transformation from true equator and equinox to Earth equator and 
//   Greenwich meridian system 
//
// Input/Output:
//
//   Mjd_UT1   Modified Julian Date UT1
//   <return>  Greenwich Hour Angle matrix
//
//  p.181
//  5.4.1 Rotation About the Celestial Ephemeris Pole
//  5.67
//------------------------------------------------------------------------------

Matrix GHAMatrix (double Mjd_UT1)
{
  return  R_z ( GAST(Mjd_UT1) );
}


//------------------------------------------------------------------------------
//
// PoleMatrix
//
// Purpose:
//
//   Transformation from pseudo Earth-fixed to Earth-fixed coordinates
//   for a given date
//
// Input/Output:
//
//   Mjd_UTC   Modified Julian Date UTC
//   <return>  Pole matrix
//
//  p.185
//  5.4.4 Transformation to the International Reference Pole
//  5.75
//------------------------------------------------------------------------------

Matrix PoleMatrix (double Mjd_UTC)
{
   return  R_y(-IERS::x_pole(Mjd_UTC)) * R_x(-IERS::y_pole(Mjd_UTC));
}


//------------------------------------------------------------------------------
//
// Geodetic (class implementation)
//
// Purpose:
//
//   Class (with all public elements) for handling geodetic coordinates
//
//  p.185-189
//  5.5 Geodetic Datums
//------------------------------------------------------------------------------

// Default constructor

Geodetic::Geodetic ()
 : lon(0.0), lat(0.0), h(0.0)
{
}

// Simple constructor

Geodetic::Geodetic (double lambda, double phi, double alt)
{
  lon=lambda; lat=phi; h=alt;
}

// Constructor for geodetic coordinates from given position

Geodetic::Geodetic (Vector r,                         // Position vector [m]
                    double R_equ,                     // Equator radius [m]
                    double f)                         // Flattening
{

  const double  eps     = 1.0e3*eps_mach;   // Convergence criterion 
  const double  epsRequ = eps*R_equ;
  const double  e2      = f*(2.0-f);        // Square of eccentricity
  
  const double  X = r(0);                   // Cartesian coordinates
  const double  Y = r(1);
  const double  Z = r(2);
  const double  rho2 = X*X + Y*Y;           // Square of distance from z-axis
  
  // Iteration 

  double  dZ, dZ_new, SinPhi;
  double  ZdZ, Nh, N;

  dZ = e2*Z;
  for(;;) {
    ZdZ    =  Z + dZ;
    Nh     =  sqrt ( rho2 + ZdZ*ZdZ ); 
    SinPhi =  ZdZ / Nh;                    // Sine of geodetic latitude
    N      =  R_equ / sqrt(1.0-e2*SinPhi*SinPhi); 
    dZ_new =  N*e2*SinPhi;
    if ( fabs(dZ-dZ_new) < epsRequ ) break;
    dZ = dZ_new;
  }
    
  // Longitude, latitude, altitude

  lon = atan2 ( Y, X );
  lat = atan2 ( ZdZ, sqrt(rho2) );
  h   = Nh - N;

}

// Position vector [m] from geodetic coordinates

Vector Geodetic::Position (double R_equ,   // Equator radius [m]
                           double f    )   // Flattening
  const
{  

  const double  e2     = f*(2.0-f);        // Square of eccentricity
  const double  CosLat = cos(lat);         // (Co)sine of geodetic latitude
  const double  SinLat = sin(lat);

  double  N;
  Vector  r(3);
      
  // Position vector 

  N = R_equ / sqrt(1.0-e2*SinLat*SinLat);

  r(0) =  (         N+h)*CosLat*cos(lon);
  r(1) =  (         N+h)*CosLat*sin(lon);
  r(2) =  ((1.0-e2)*N+h)*SinLat;

  return r;

}

// Transformation to local tangent coordinates

Matrix Geodetic::LTC_Matrix () const
{
  return LTCMatrix(lon,lat);
}


//------------------------------------------------------------------------------
//
// LTCMatrix
//
// Purpose:
//
//   Transformation from Greenwich meridian system to local tangent coordinates
//
// Input/Output:
//
//   lambda    Geodetic East longitude [rad]
//   phi       Geodetic latitude [rad]
//   <return>  Rotation matrix from the Earth equator and Greenwich meridian
//             to the local tangent (East-North-Zenith) coordinate system
//
//------------------------------------------------------------------------------

Matrix LTCMatrix (double lambda, double phi)
{
  
  Matrix  M(3,3);
  double  Aux;
  
  // Transformation to Zenith-East-North System
  M = R_y(-phi)*R_z(lambda);
  
  // Cyclic shift of rows 0,1,2 to 1,2,0 to obtain East-North-Zenith system
  for (int j=0; j<3; j++) {
    Aux=M(0,j); M(0,j)=M(1,j); M(1,j)=M(2,j); M(2,j)= Aux;
  }
  
  return  M;
  
}


//------------------------------------------------------------------------------
//
// AzEl
//
// Purpose:
//
//   Computes azimuth and elevation from local tangent coordinates
//
// Input/Output:
//
//   s   Topocentric local tangent coordinates (East-North-Zenith frame)
//   A   Azimuth [rad]
//   E   Elevation [rad]
//
//------------------------------------------------------------------------------

void AzEl ( const Vector& s, double& A, double& E ) 
{
  A = atan2(s(0),s(1));
  A = ((A<0.0)? A+pi2 : A);
  E = atan ( s(2) / sqrt(s(0)*s(0)+s(1)*s(1)) );
}


//------------------------------------------------------------------------------
//
// AzEl
//
// Purpose:
//
//   Computes azimuth, elevation and partials from local tangent coordinates
//
// Input/Output:
//
//   s      Topocentric local tangent coordinates (East-North-Zenith frame)
//   A      Azimuth [rad]
//   E      Elevation [rad]
//   dAds   Partials of azimuth w.r.t. s
//   dEds   Partials of elevation w.r.t. s
//
//------------------------------------------------------------------------------

void AzEl( const Vector& s, double& A, double& E, Vector& dAds, Vector& dEds ) 
{
  const double rho = sqrt(s(0)*s(0)+s(1)*s(1));
  // Angles
  A = atan2(s(0),s(1));
  A = ((A<0.0)? A+pi2 : A);
  E = atan ( s(2) / rho );
  // Partials
  dAds = Vector( s(1)/(rho*rho), -s(0)/(rho*rho), 0.0 );
  dEds = Vector( -s(0)*s(2)/rho, -s(1)*s(2)/rho , rho ) / Dot(s,s);
}

Die Fehler-Ausgabe sieht so aus:
Code:
SAT_RefSys1.cpp: In function ‘Matrix EclMatrix(double)’:
SAT_RefSys1.cpp:172: error: ‘R_x’ was not declared in this scope
SAT_RefSys1.cpp: In function ‘Matrix PrecMatrix(double, double)’:
SAT_RefSys1.cpp:217: error: ‘R_z’ was not declared in this scope
SAT_RefSys1.cpp:217: error: ‘R_y’ was not declared in this scope
SAT_RefSys1.cpp: In function ‘Matrix NutMatrix(double)’:
SAT_RefSys1.cpp:437: error: ‘R_x’ was not declared in this scope
SAT_RefSys1.cpp:437: error: ‘R_z’ was not declared in this scope
SAT_RefSys1.cpp: In function ‘Matrix NutMatrixSimple(double)’:
SAT_RefSys1.cpp:487: error: ‘R_x’ was not declared in this scope
SAT_RefSys1.cpp:487: error: ‘R_z’ was not declared in this scope
SAT_RefSys1.cpp: In function ‘Matrix GHAMatrix(double)’:
SAT_RefSys1.cpp:620: error: ‘R_z’ was not declared in this scope
SAT_RefSys1.cpp: In function ‘Matrix PoleMatrix(double)’:
SAT_RefSys1.cpp:645: error: ‘R_y’ was not declared in this scope
SAT_RefSys1.cpp:645: error: ‘R_x’ was not declared in this scope
SAT_RefSys1.cpp: In function ‘Matrix LTCMatrix(double, double)’:
SAT_RefSys1.cpp:773: error: ‘R_y’ was not declared in this scope
SAT_RefSys1.cpp:773: error: ‘R_z’ was not declared in this scope
make[4]: *** [SAT_RefSys1.o] Error 1

Ich hoffe nicht, dass es zu weiteren Fehlern in den (anderen) Quellen kommt, da die urspr. Programmierer dieses Codes nicht mehr erreichbar sind.

MfG,
upperlimit
 
Wenn Du das Programm einfach nur kompilieren willst, ist das einfachste sicherlich, wenn Du Dir gcc-3.3.5 runterlaedst (gcc.gnu.org), kompilierst und damit dann das Programm kompilierst.

Wenn Dir an dem Code was gelegen ist, musst Du die Fehlermeldungen Stueck fuer Stueck durcharbeiten und beheben. Suche mal nach der Funktion R_x und wo sie deklariert ist. Ich habe sie in den Dateien, die Du geschickt hast, nicht gefunden. Sie muss an der Stelle, an der sie in SAT_RefSys1.cpp zum ersten Mal benutzt wird, bereits bekannt, d.h., deklariert sein. In der SAT_VecMat.h ist sie als friend zur Klasse Matrix eingetragen, doch ist das keine Deklaration, soweit ich weiss.
 
rikola schrieb:
Ich wuerde mal sagen, dass der Code nichts Versionsspezifisches enthaelt. gcc-4.0 ist einfach strenger und enger am Standard. In der Datei SAT_VecMat.h wird die Klasse Matrix innerhalb der Klasse Vector benutzt, ohne sie vorher zu deklarieren. Schreibe einfach ueber Zeile 44, also vor der Deklaration der Klasse Vector, ein 'class Matrix;' hin, dann sollte zumindest diese Fehlermeldung verschwinden. Es wundert mich, dass gcc-3.3 das hat durchgehen lassen.
Och ich kenne einige für gcc 3.3.X geschriebene Programme, die sich nicht so einfach mit gcc 4 compilieren lassen.

Da sind eben Anpassungen nötig.
Gruß Wolfgang
 
Wolfgang schrieb:
Och ich kenne einige für gcc 3.3.X geschriebene Programme, die sich nicht so einfach mit gcc 4 compilieren lassen.

Da sind eben Anpassungen nötig.
Gruß Wolfgang
Schon klar. Ich bezog mich auf den von upperlimt angegebenen Code und v.a. die Fehlermeldung, die einfach typisch fuer eine fehlende Deklaration ist.
 
Hallo,

inzwischen habe ich den Code mit gcc 3.4.6 erfolgreich getestet, negativ mit gcc 4.1.1.


Sofern rikola Recht hat, sollte es doch moeglich sein, die strengen ISO-Vorschriften vom gcc 4.1.0 mittels eines Schalters/Option zu entkraeften, so dass auch Code fuer aeltere gcc-Versionen problemlos kompiliert werden kann, richtig ?
Ich habe bisher '-std=c++98', '-std=gnu++98' sowie '-ansi' ausprobiert, erfolglos.

@rikola,
Du hast recht, die Deklaration von R_x, R_y und R_z erfolgt in einer anderen Datei, genannt SAT_VecMat1.cpp (s. Code)

Code:
//------------------------------------------------------------------------------
//
// SAT_VecMat1.cpp
// 
// Purpose: 
//
//   Vector/matrix operations
//   
// Notes:
//
//   This software is protected by national and international copyright. 
//   Any unauthorized use, reproduction or modificaton is unlawful and 
//   will be prosecuted. Commercial and non-private application of the 
//   software in any form is strictly prohibited unless otherwise granted
//   by the authors.
//
//   The code is provided without any warranty; without even the implied 
//   warranty of merchantibility or fitness for a particular purpose.
//
// Last modified:
//
//   2000/03/04  OMO  Final version (1st edition)
//
// (c) 1999-2000  O. Montenbruck, E. Gill
//
//------------------------------------------------------------------------------

#include <cmath>
#include <iostream>
#include <iomanip>

#include "SAT_VecMat.h"

/*#ifdef __GNUC__  // GNU C++ adaptation
#include "GNU_iomanip.h"
#endif*/

using namespace std;


//------------------------------------------------------------------------------
// 
// LU_Decomp
//
// Purpose:
//
//   LU-Decomposition.
//
//   Given an nxn matrix A, this routine replaces it by the LU decomposition
//   of a rowwise permutation of itself. A is output, arranged as in 
//   equation (2.3.14) of Press et al. (1986); Indx is an ouput vector which
//   records the row permutation effected by partial pivoting. This routine is 
//   used in combination with LU_BackSub to solve linear equations or invert 
//   a matrix.
//
// Input/output:
//
//   A       Square matrix; replaced by LU decomposition of permutation of A
//           on output
//   Indx    Permutation index vector
//
// Note:
//
//   Adapted from LUDCMP of Press et al. (1986).
//
//------------------------------------------------------------------------------

void LU_Decomp ( Matrix& A, Vector& Indx )
{

  // Constants
      
  const int    n    = A.size1();
  const double tiny = 1.0e-20;       // A small number

  // Variables
      
  int     i,j,imax,k;
  double  aAmax, Sum, Dum;
  Vector  V(n); 

  // Loop over rows to get scaling information

  for (i=0; i<n; i++) {
    aAmax = 0.0;
    for (j=0;j<n;j++) if (fabs(A(i,j)) > aAmax ) aAmax=fabs(A(i,j));
    if (aAmax==0.0) {
      // No nonzero largest element
      cerr << "ERROR: Singular matrix A in LU_Decomp";
      exit(1);
    };
    V(i) = 1.0/aAmax;           // V stores the implicit scaling of each row
  };

  // Loop over columns of Crout's method
      
  for ( j=0; j<n; j++ ) {
      
    if (j > 0) {
      for ( i=0; i<j; i++ ) {   // This is equation 2.3.12 except for i=j
        Sum = A(i,j);
        if (i>0) {
          for ( k=0; k<i; k++ )  Sum -= A(i,k)*A(k,j);
          A(i,j) = Sum;
        };
      };
    };
      
    aAmax=0.0;                  // Initialize for the search of the largest
                                // pivot element
      
    for ( i=j; i<n; i++ ) {     // This is i=j of equation 2.3.12 and 
      Sum = A(i,j);             // i=j+1..N of equation 2.3.13
      if (j > 0) {
        for ( k=0; k<j; k++ ) Sum -= A(i,k)*A(k,j);
        A(i,j) = Sum;
      };
      Dum = V(i)*fabs(Sum);     // Figure of merit for the pivot
      if (Dum >= aAmax) {       // Is it better than the best so far ?
        imax  = i;
        aAmax = Dum;
      };
    };
      
    if (j != imax) {            // Do we need to interchange rows?
      for ( k=0; k<n; k++) {    // Yes, do so ...
        Dum = A(imax,k);
        A(imax,k) = A(j,k);
        A(j,k) = Dum;
      }
      V(imax) = V(j);           // Also interchange the scale factor 
    };
      
    Indx(j) = imax;
      
    if (j != n-1) {             // Now finally devide by the pivot element
      if (A(j,j) == 0.0) {      // If the pivot element is zero the matrix 
        A(j,j) = tiny;          // is singular (at least to the precision of
      };                        // the algorithm). For some applications on
      Dum=1.0/A(j,j);           // singular matrices, it is desirable to 
      for (i=j+1;i<n;i++) {     // substitude tiny for zero. 
        A(i,j)=A(i,j)*Dum;
      };
    };

  };   // Go back for the next column in the reduction

  if (A(n-1,n-1)==0.0) A(n-1,n-1)=tiny; 

};


//------------------------------------------------------------------------------
//
// LU_BackSub
//
// Purpose:
//
//   LU Backsubstitution
//
//   Solves the set of n linear equations Ax=b. Here A is input, not as the 
//   matrix A but rather as its LU decomposition, determined by the function
//   LU_Decomp. b is input as the right-hand side vector b, and returns with
//   the solution vector x. A and Indx are not modified by this function and 
//   can be left in place for successive calls with different right-hand 
//   sides b. This routine takes into account the posssibility that B will  
//   begin with many zero elements, so it is efficient for use in matrix
//   inversions.
//
// Input/output:
//
//   A       LU decomposition of permutation of A
//   Indx    Permutation index vector
//   b       Right-hand side vector b; replaced by solution x of Ax=b on output
//
//------------------------------------------------------------------------------

void LU_BackSub ( Matrix& A, Vector& Indx, Vector& b )
{ 

  // Constants
      
  const int  n = A.size1();

  // Local variables

  int     ii,i,ll,j;
  double  Sum;

  //
  // Start
  //

  ii = -1;                      // When ii is set to a nonegative value, it will
                                // become the first nonvanishing element of B. 
  for (i=0; i<n; i++) {         // We now do the forward substitution.
    ll = (int) Indx(i);         // The only wrinkle is to unscramble the 
    Sum = b(ll);                // permutation as we go.
    b(ll) = b(i);
    if (ii != -1) {
      for (j=ii; j<i; j++) Sum -= A(i,j)*b(j);
    }
    else {
      if (Sum != 0.0) ii = i;   // A nonzero element was encountered, so from 
    };                          // now on we will have to do the sums in the
    b(i) = Sum;                 // loop above.
   };
      
   for (i=n-1; i>=0; i--) {     // Now we do the backsubstitution, eqn 2.3.7.
     Sum=b(i);
     if (i<n-1) {
       for (j=i+1;j<n;j++) {
          Sum = Sum-A(i,j)*b(j);
       };
     };
     b(i) = Sum/A(i,i);         // Store a component of the solution vector X.
   };

};



//------------------------------------------------------------------------------
//
// Vector (class implementation)
//
// Purpose:
//
//   Vector data type and associated operations
//
//------------------------------------------------------------------------------


// Constructors, destructor

Vector::Vector ()                              // Vector without elements
  : n(0) 
{
  v = 0;
}

Vector::Vector (int Size)                      // Creates null vector
 : n(Size)
{
  v = new double [Size];
  for (int i=0; i<Size; i++) v[i]=0.0;
}

Vector::Vector (const Vector& V)               // Vector copy
 : n(V.n)
{
  v = new double [V.n];
  for (int i=0; i<V.n; i++) v[i]=V.v[i];
}

Vector::Vector (const double* p, int N)        // Array copy
  : n(N)
{
  v = new double [N];
  for (int i=0; i<N; i++) v[i]=p[i];
}

Vector::Vector (double x, double y, double z)  // 3dim-Vector
  : n(3)
{
  v = new double [3];
  v[0]=x; v[1]=y; v[2]=z;
}

Vector::Vector (double x, double y, double z,   // 6dim-Vector
                double X, double Y, double Z)  
  : n(6)
{
  v = new double [6];
  v[0]=x; v[1]=y; v[2]=z;
  v[3]=X; v[4]=Y; v[5]=Z;
}

Vector::~Vector() 
{ 
  delete [] v; 
}

// Component access

Vector Vector::slice (int first, int last) const
{
  Vector Aux(last-first+1);
  for (int i=first; i<=last; i++) Aux.v[i-first]=v[i];
  return Aux;
}


// Square root of vector elements

Vector Vector::Sqrt()
{
  Vector Aux(n);
  for (int i=0; i<n; i++) Aux.v[i]=sqrt(v[i]);
  return Aux;
}


// Assignment

Vector& Vector::operator=(const double value)
{
  for (int i=0; i<n; i++) v[i]=value;
  return (*this);
}

Vector& Vector::operator=(const Vector& V)
{
  if (this == &V) return (*this);
  // Allocate vector if still empty
  if (v==0) {
    n = V.n; 
    v = new double [V.n];
  };
  // Check dimension
  if (n!=V.n) {
    cerr << "ERROR: Incompatible sizes in Vector operator=(Vector)" << endl;
    exit(1);
  };
  // Copy elements
  for (int i=0; i<n; i++) v[i]=V.v[i];
  return (*this);
}


// Concatenation

Vector Stack (const Vector& a, const Vector& b)
{
  int    i;
  Vector c(a.size()+b.size());
  for (i=0;i<a.size();i++) c(i)=a(i);
  for (i=0;i<b.size();i++) c(i+a.size())=b(i);
  return c;
}


// Vector from polar angles

Vector VecPolar (double azim, double elev, double r)
{
  return Vector(r*cos(azim)*cos(elev),r*sin(azim)*cos(elev),r*sin(elev));
}


// Vector addition/subtraction with assignment

void Vector::operator += (const Vector& V)
{
  if (n!=V.n) {
    cerr << "ERROR: Incompatible shape in Vector operator+=(Vector)" << endl;
    exit(1);
  };
  for (int i=0; i<n; i++) v[i]+=V.v[i];
}

void Vector::operator -= (const Vector& V)
{
  if (n!=V.n) {
    cerr << "ERROR: Incompatible shape in Vector operator-=(Vector)" << endl;
    exit(1);
  };
  for (int i=0; i<n; i++) v[i]-=V.v[i];
}


// Dot product, norm, cross product

double Dot (const Vector& left, const Vector& right)
{
  if (left.n!=right.n) {
    cerr << "ERROR: Incompatible shape in Dot(Vector,Vector)" << endl;
    exit(1);
  };
  double Sum = 0.0;
  for (int i=0; i<left.n; i++) Sum+=left.v[i]*right.v[i];
  return Sum;
}

double Norm (const Vector& V)
{
  return sqrt(Dot(V,V));
}

Vector Cross (const Vector& left, const Vector& right)
{
  if ( (left.n!=3) || (right.n!=3) ) {
    cerr << "ERROR: Invalid dimension in Cross(Vector,Vector)" << endl;
    exit(1);
  };
  Vector Result(3);
  Result.v[0] = left.v[1]*right.v[2] - left.v[2]*right.v[1];
  Result.v[1] = left.v[2]*right.v[0] - left.v[0]*right.v[2];
  Result.v[2] = left.v[0]*right.v[1] - left.v[1]*right.v[0];
  return Result;
}


// Scalar multiplication and division of a vector

Vector operator * (double value, const Vector& V)
{
  Vector Aux(V.n);
  for (int i=0; i<V.n; i++) Aux.v[i]=value*V.v[i];
  return Aux;
}

Vector operator * (const Vector& V, double value)
{
  return value*V;
}

Vector operator / (const Vector& V, double value)
{
  Vector Aux(V.n);
  for (int i=0; i<V.n; i++) Aux.v[i]=V.v[i]/value;
  return Aux;
}


// Negation of a vector (unary minus)

Vector operator - (const Vector& V)
{
  Vector Aux(V.n);
  for (int i=0; i<V.n; i++) Aux.v[i]=-V.v[i];
  return Aux;
}


// Vector addition and subtraction

Vector operator + (const Vector& left, const Vector& right)
{
  if (left.n!=right.n) {
    cerr << "ERROR: Incompatible shape in +(Vector,Vector)" << endl;
    exit(1);
  };
  Vector Aux(left.n);
  for (int i=0; i<left.n; i++) Aux.v[i]=left.v[i]+right.v[i];
  return Aux;
}

Vector operator - (const Vector& left, const Vector& right)
{
  if (left.n!=right.n) {
    cerr << "ERROR: Incompatible shape in -(Vector,Vector)" << endl;
    exit(1);
  };
  Vector Aux(left.n);
  for (int i=0; i<left.n; i++) Aux.v[i]=left.v[i]-right.v[i];
  return Aux;
}

// Vector output

ostream& operator << (ostream& os, const Vector& Vec)
{
  int w = os.width();
  for (int i=0; i<Vec.size(); i++)
    os << setw(w) << Vec(i);
//os << endl;
  return os;
}



//------------------------------------------------------------------------------
//
// Matrix (class implementation)
//
// Purpose:
//
//   Matrix data type and associated operations
//
//------------------------------------------------------------------------------


// Constructors

Matrix::Matrix ()                          // Matrix without elements
  : n(0), m(0) 
{
  M = 0;
}

Matrix::Matrix (int dim1, int dim2)        // Nullmatrix of specified shape 
  : n(dim1), m(dim2)
{
  int i,j;
  // Memory allocation
  M = new double*[dim1];
  for (i=0; i<dim1; i++) M[i] = new double[dim2];
  // Initialization
  for (i=0; i<dim1; i++) {
    for (j=0; j<dim2; j++) M[i][j]=0.0;
  }
}

Matrix::Matrix (const Matrix& M_)          // Copy
{
  int i,j;
  n = M_.n;
  m = M_.m;
  // Memory allocation
  M = new double*[n];
  for (i=0; i<n; i++) M[i] = new double[m];
  // Initialization
  for (i=0; i<n; i++) {
    for (j=0; j<m; j++) M[i][j]=M_.M[i][j];
  }
}

Matrix::Matrix (const double* p, int dim1, int dim2)   // Array copy
{
  int i,j;
  n = dim1;
  m = dim2;
  // Memory allocation
  M = new double*[n];
  for (i=0; i<n; i++) M[i] = new double[m];
  // Initialization
  for (i=0; i<n; i++) {
    for (j=0; j<m; j++) M[i][j]=p[i*dim2+j];
  }
  
}

Matrix::~Matrix() 
{ 
  for (int i=0; i<n; i++) delete[] M[i];
  delete [] M; 
};


// Assignment

Matrix& Matrix::operator=(const double value)
{
  for (int i=0; i<n; i++) 
    for (int j=0; j<n; j++) 
      M[i][j]=value;
  return (*this);
}

Matrix& Matrix::operator=(const Matrix& M_)
{
  int i,j;
  if (this == &M_) return (*this);
  // Allocate matrix if still empty
  if (M==0) {
    n = M_.n; 
    m = M_.m;
    M = new double*[n];
    for (i=0; i<n; i++) M[i] = new double[m];
  };
  if ( (n!=M_.n) || (m!=M_.m) ) {
    cerr << "ERROR: Incompatible shapes in Matrix operator=(Matrix)" << endl;
    exit(1);
  };
  for (i=0; i<n; i++) 
    for (j=0; j<m; j++) 
      M[i][j]=M_.M[i][j];
  return (*this);
}


// Component access

Vector Matrix::Col(int j) const
{
  Vector Res(n);
  for (int i=0; i<n; i++)  Res.v[i]=M[i][j];
  return Res;
}

Vector Matrix::Row(int i) const
{
  Vector Res(m);
  for (int j=0; j<m; j++)  Res.v[j]=M[i][j];
  return Res;
}

Vector Matrix::Diag() const
{
  if (n!=m) {
    cerr << "ERROR: Invalid shape in Matrix.Diag()" << endl;
    exit(1);
  };
  Vector Vec(n);
  for (int i=0; i<n; i++) Vec.v[i] = M[i][i];
  return Vec;
}


Matrix Matrix::slice(int first_row, int last_row, int first_col, int last_col)
{
  if (first_row<0 || last_row<first_row || n-1<last_row ||
      first_col<0 || last_col<first_col || m-1<last_col) {
    cerr << "ERROR: Invalid arguments in Matrix.slice()" << endl;
    exit(1);
  };
  Matrix Aux(last_row-first_row+1,last_col-first_col+1);
  for (int i=0;i<=last_row-first_row;i++)
    for (int j=0;j<=last_col-first_col;j++)
       Aux(i,j) = M[i+first_row][j+first_col];
  return Aux;
}


void Matrix::SetCol(int j, const Vector& Col) 
{
  if (Col.size()!=n) {
    cerr << "ERROR: Incompatible shapes in Matrix.SetCol()" << endl;
    exit(1);
  };
  if (j<0 || m<=j) {
    cerr << "ERROR: Column index out of range in Matrix.SetCol()" << endl;
    exit(1);
  };
  for (int i=0; i<n; i++) M[i][j]=Col(i);
}

void Matrix::SetRow(int i, const Vector& Row)
{
  if (Row.size()!=m) {
    cerr << "ERROR: Incompatible shapes in Matrix.SetRow()" << endl;
    exit(1);
  };
  if (i<0 || n<=i) {
    cerr << "ERROR: Row index out of range in Matrix.SetRow()" << endl;
    exit(1);
  };
  for (int j=0; j<n; j++) M[i][j]=Row(j);
}


// Unit matrix

Matrix Id(int Size)
{
  Matrix Aux(Size,Size);      
  for (int i=0; i<Size; i++) Aux.M[i][i] = 1.0;
  return Aux;
}


// Diagonal matrix

Matrix Diag(const Vector& Vec)
{
  Matrix Mat(Vec.n,Vec.n);
  for (int i=0; i<Vec.n; i++) Mat.M[i][i] = Vec.v[i];
  return Mat;
}


// Elementary rotations

Matrix R_x(double Angle)
{
  const double C = cos(Angle);
  const double S = sin(Angle);
  Matrix U(3,3);
  U.M[0][0] = 1.0;  U.M[0][1] = 0.0;  U.M[0][2] = 0.0;
  U.M[1][0] = 0.0;  U.M[1][1] =  +C;  U.M[1][2] =  +S;
  U.M[2][0] = 0.0;  U.M[2][1] =  -S;  U.M[2][2] =  +C;
  return U;
}

Matrix R_y(double Angle)
{
  const double C = cos(Angle);
  const double S = sin(Angle);
  Matrix U(3,3);
  U.M[0][0] =  +C;  U.M[0][1] = 0.0;  U.M[0][2] =  -S;
  U.M[1][0] = 0.0;  U.M[1][1] = 1.0;  U.M[1][2] = 0.0;
  U.M[2][0] =  +S;  U.M[2][1] = 0.0;  U.M[2][2] =  +C;
  return U;
}

Matrix R_z(double Angle)
{
  const double C = cos(Angle);
  const double S = sin(Angle);
  Matrix U(3,3);
  U.M[0][0] =  +C;  U.M[0][1] =  +S;  U.M[0][2] = 0.0;
  U.M[1][0] =  -S;  U.M[1][1] =  +C;  U.M[1][2] = 0.0;
  U.M[2][0] = 0.0;  U.M[2][1] = 0.0;  U.M[2][2] = 1.0;
  return U;
}


// Transposition

Matrix Transp(const Matrix& Mat)
{
  Matrix T(Mat.m,Mat.n);
  for ( int i=0; i<T.n; i++ )
    for ( int j=0; j<T.m; j++ )
      T.M[i][j] = Mat.M[j][i];
  return T;
}


// Inverse

Matrix Inv(const Matrix& Mat)
{
  const int n = Mat.n;

  Matrix LU(n,n), Inverse(n,n);
  Vector b(n), Indx(n);

  if (Mat.m!=Mat.n) {
    cerr << "ERROR: Invalid shape in Inv(Matrix)" << endl;
    exit(1);
  };

  // LU decomposition 

  LU = Mat;
  LU_Decomp ( LU, Indx );

  // Solve Ax=b for  unit vectors b_1..b_n

  for (int j=0; j<n; j++ ) {
    b=0.0; b(j)= 1.0;                     // Set b to j-th unit vector 
    LU_BackSub ( LU, Indx, b );           // Solve Ax=b 
    Inverse.SetCol(j,b);                  // Copy result
  };

  return Inverse;

}


// Scalar multiplication and division of a vector

Matrix operator * (double value, const Matrix& Mat)
{
  Matrix Aux(Mat.n,Mat.m);
  for (int i=0; i<Mat.n; i++) 
    for (int j=0; j<Mat.m; j++) 
      Aux.M[i][j]=value*Mat.M[i][j];
  return Aux;
}

Matrix operator * (const Matrix& Mat, double value)
{
  return value*Mat;
}

Matrix operator / (const Matrix& Mat, double value)
{
  Matrix Aux(Mat.n,Mat.m);
  for (int i=0; i<Mat.n; i++) 
    for (int j=0; j<Mat.m; j++) 
      Aux.M[i][j]=Mat.M[i][j]/value;
  return Aux;
}


// Unary minus

Matrix operator - (const Matrix& Mat)
{
  Matrix Aux(Mat.n,Mat.m);
  for (int i=0; i<Mat.n; i++) 
    for (int j=0; j<Mat.m; j++) 
      Aux.M[i][j]=-Mat.M[i][j];
  return Aux;
}


// Matrix addition and subtraction

Matrix operator + (const Matrix& left, const Matrix& right)
{
  if ( (left.n!=right.n) || (left.m!=right.m) ) {
    cerr << "ERROR: Incompatible shape in +(Matrix,Matrix)" << endl;
    exit(1);
  };
  Matrix Aux(left.n,left.m);
  for (int i=0; i<left.n; i++) 
    for (int j=0; j<left.m; j++) 
      Aux.M[i][j] = left.M[i][j] + right.M[i][j];
  return Aux;
}

Matrix operator - (const Matrix& left, const Matrix& right)    
{
  if ( (left.n!=right.n) || (left.m!=right.m) ) {
    cerr << "ERROR: Incompatible shape in -(Matrix,Matrix)" << endl;
    exit(1);
  };
  Matrix Aux(left.n,left.m);
  for (int i=0; i<left.n; i++) 
    for (int j=0; j<left.m; j++) 
      Aux.M[i][j] = left.M[i][j] - right.M[i][j];
  return Aux;
}


// Matrix product

Matrix operator * (const Matrix& left, const Matrix& right)
{
  if (left.m!=right.n) {
    cerr << "ERROR: Incompatible shape in *(Matrix,Matrix)" << endl;
    exit(1);
  };
  Matrix Aux(left.n,right.m);
  double Sum;
  for (int i=0; i<left.n; i++) 
    for (int j=0; j<right.m; j++) {
      Sum = 0.0;
      for (int k=0; k<left.m; k++) 
        Sum += left.M[i][k] * right.M[k][j];
      Aux.M[i][j] = Sum;
    }
  return Aux;
}


// Vector/matrix product

Vector operator * (const Matrix& Mat, const Vector& Vec)
{
  if (Mat.m!=Vec.n) {
    cerr << "ERROR: Incompatible shape in *(Matrix,Vector)" << endl;
    exit(1);
  };
  Vector Aux(Mat.n);
  double Sum;
  for (int i=0; i<Mat.n; i++) {
    Sum = 0.0;
    for (int j=0; j<Mat.m; j++) 
      Sum += Mat.M[i][j] * Vec.v[j];
    Aux.v[i] = Sum;
  }
  return Aux;
}

Vector operator * (const Vector& Vec, const Matrix& Mat)
{
  if (Mat.n!=Vec.n) {
    cerr << "ERROR: Incompatible shape in *(Vector,Matrix)" << endl;
    exit(1);
  };
  Vector Aux(Mat.m);
  double Sum;
  for (int j=0; j<Mat.m; j++) {
    Sum = 0.0;
    for (int i=0; i<Mat.n; i++) 
      Sum += Vec.v[i] * Mat.M[i][j];
    Aux.v[j] = Sum;
  }
  return Aux;
}


// Dyadic product

Matrix Dyadic (const Vector& left, const Vector& right)
{
  Matrix Mat(left.n,right.n);
  for (int i=0;i<left.n;i++)
    for (int j=0;j<right.n;j++)
      Mat.M[i][j] = left.v[i]*right.v[j];
  return Mat;
}


// Matrix output

ostream& operator << (ostream& os, const Matrix& Mat)
{
  int w = os.width();
  for (int i=0; i<Mat.size1(); i++) {
    for (int j=0; j<Mat.size2(); j++)
      os << setw(w) << Mat(i,j);
    os << endl;
  }
  return os;
}


MfG,
upperlimit
 
naja mit der option -std=c++98 kannst du ja auch nicht viel erreichen.
Damit sagst du ja du willst den c++ standard von 98, dies ist der aktuelle
ISO-standard. Den willst du ja grade nicht!
um ehrlich zu sein, kann ich dir leider auch keinen schalter nennen, der den
standard "abschlatet", btw denke ich auch, dass dies vll nicht wirklich der
richtige weg ist.
 

Ähnliche Themen

Shellskript - Fehler in Cron

Aufgabe in C

Amarok stoppt nach jedem Stück

Samba 4 Gast Zugang unter Ubuntu funktioniert nicht

Problem mit HSPA+ Modem Huawei E353 - Installation unmöglich?

Zurück
Oben