00001
00006 #include "bpm/bpm_nr.h"
00007
00008
00009
00010 complex_t complex( double re, double im ) {
00011 complex_t z;
00012 z.re = re;
00013 z.im = im;
00014 return z;
00015 }
00016
00017 double c_real( complex_t z ) { return z.re; }
00018 double c_imag( complex_t z ) { return z.im; }
00019 double c_abs( complex_t z ) { return sqrt( z.re*z.re+z.im*z.im ); }
00020 double c_arg( complex_t z ) { return atan2( z.im, z.re ); }
00021
00022 complex_t c_conj( complex_t z ) { return complex( z.re, -z.im ); }
00023 complex_t c_neg( complex_t z ) { return complex( -z.re, -z.im ); }
00024
00025 complex_t c_sum( complex_t z1, complex_t z2 ) { return complex( z1.re + z2.re, z1.im + z2.im ); }
00026 complex_t c_diff( complex_t z1, complex_t z2 ) { return complex( z1.re - z2.re, z1.im - z2.im ); }
00027
00028 complex_t c_mult( complex_t z1, complex_t z2 ) {
00029 return complex( z1.re*z2.re - z1.im*z2.im,
00030 z1.re*z2.im + z1.im*z2.re );
00031 }
00032
00033 complex_t c_scale( double r, complex_t z ) { return complex( r*z.re, r*z.im); }
00034
00035
00036 complex_t c_div( complex_t z1, complex_t z2 ) {
00037 double r = c_norm2( z2 );
00038 return complex( ((z1.re*z2.re) + (z1.im*z2.im)) / r,
00039 ((z1.im*z2.re) - (z1.re*z2.im)) / r );
00040 }
00041
00042 complex_t c_sqr( complex_t z ) { return c_mult( z, z ); }
00043
00044 double c_norm2( complex_t z ) { return ( z.re*z.re + z.im*z.im ); }
00045
00046 complex_t c_exp( complex_t z ) { return complex( exp(z.re)*cos(z.im), exp(z.re)*sin(z.im) ); }
00047
00048 complex_t c_sqrt( complex_t z ) {
00049 double r = c_abs( z );
00050 complex_t u = complex( sqrt( (r + z.re ) / 2. ), sqrt( ( r - z.re ) / 2. ) );
00051 if( z.im < 0.0 ) u.im = -u.im;
00052 return u;
00053 }
00054
00055
00056 int c_isequal( complex_t z1, complex_t z2 ) {
00057 if ( ( z1.re == z2.re ) && ( z1.im == z2.im ) ) {
00058 return 1;
00059 } else return 0;
00060 }
00061