jon@hanauma (Jon Claerbout) (06/25/89)
Scientific programming in C++ ?
ABSTRACT:
I converted some simple scientific Fortran/Ratfor programs into C++
to see if they would look suitable for a textbook
such as my last book "Imaging the Earth's Interior".
I conclude C++ is about as good as ratfor.
Unfortunately, mixing Fortran with C++
ranges from undocumented to impossible.
BOOK REPORT
The new Lippman C++ book looks like a replacement for the Stroustrup book
since it fully describes the new AT&T version 2.0.
Pedagogically it is a big improvement too.
Since both Stroustrup and Lippman describe both C and C++
I infer they mean eventually to replace C by C++
(else why the 100+ extra pages to explain C
which K&R already do beautifully)?
A section called "linkage to other languages" mentions C but not Fortran.
PROGRAMS
We couldn't link gnu C++ mains to fortran subroutines.
I converted two simple fortran scientific programs to a C++ style
designed to please fortran users. Here they are:
-----------------
/*
* Demonstrate a scientific C++ program in a fortran-like style
* such as would be suitable for a numerical analysis library
* or an engineering textbook with tutorial programs.
* Works on gnu c++
*/
void filter( // convolution routine.
const int nb , // number of coefficients in filter
const float bb[] , // numerator moving average filter
const int npq , // length of time series, both of pp and qq
const float pp[] , // input time series
float qq[] ) // output time series
{
int ib, i ;
for( i=0 ; i < npq ; i++ ) { // This means "do i=0,npq-1"
qq[i] = 0. ;
}
for( ib=0 ; ib < nb; ib++ ) {
for( i=0 ; i < npq ; i++ ) {
qq[i+ib] = qq[i+ib] + pp[i] * bb[ib] ;
}
}
}
// Test program follows below.
#include <stream.h>
const npq = 4;
const nb = 2;
float bb[ nb] = { 1., 1.};
float pp[npq] = { 1., 0., 1., 1. };
float qq[npq];
main() {
// void filter( int, float[], int, float[], float[]);
filter( nb, bb, npq, pp, qq);
cout << "Should print 1 1 1 2 next.\n";
for( int i=0; i<npq; i++ )
cout << qq[i] << "\n";
return 0; // return "no error"
}
----------------------------------------- another example
// Illustrate C++ with complex arithmetic.
// (Compute spectrum via slow Fourier transform.)
#include <CC/complex.h>
void spectrum(
const int nb , // length of time series
const float bb[] , // time series
const int ns , // number of points in the amplitude spectrum
float ss[] ) // amplitude spectrum
{
int ib, is ;
float top;
complex cz, cb ;
for( is=0 ; is < ns ; is++ ) { // this is "do is=0,ns-1"
cz = exp( complex( 0., (3.14159256/(ns-1)) * (is) ) ) ;
cb = bb[nb-1] ;
for( ib=nb-1 ; ib > 0 ; ib-- ) {
cb = cb * cz + bb[ib-1] ;
}
top = real( cb * conj(cb)) ;
ss[is] = sqrt( top ) ;
}
}