Download Problem Set 2 with Solutions - Computational Physics | ASTR 415 and more Assignments Astronomy in PDF only on Docsity!
ASTR415 Spring 2009 Problem Set
Solution
- My code listing is given at the end, with output for single and double precision (the macro FLOAT is used to choose precision). In the single precision case, the smaller root as computed from the quadratic formula differs from that computed from the larger root because the discriminant b^2 − 4 ac nearly equals b^2. The resulting subtraction in the formula introduces round-off error at the 20% level. With double precision the round-off error is greatly reduced.
- My code listing is given at the end, with output for single and double precision (the macro FLOAT is used to choose precision). The round-off error in the recursion formula is not random: it oscillates with steadily increasing amplitude above and below the correct solution. With double precision the round-off error is greatly reduced in the specified range for n.
- My code listing is given at the end. The code reads in the specified data file in double precision until EOF is reached (or the arrays are full). The center-of-mass position and velocity is subtracted from each particle and the total angular momentum and inertia tensor are computed. The single-precision NRiC function gaussj() is used to solve for the spin vector, which is then converted into a spin period in hours.
(a) Analysis of the supplied file is given at the end. The spin period is 4.30458 h. (b) Views projected along the Cartesian axes are given at the end, along with the sm script used to generate them.
Code Listings & Output
** roots.c */
#include #include
#define FLOAT float /* could also do "typedef float FLOAT;" */
#define SQ(x) ((x)(x)) / handy macro */
int main(void) { FLOAT x,x1,x2,xmax,xmin;
(void) printf("%s-precision version\n", sizeof(FLOAT)==sizeof(float)?"single":"double");
x = sqrt(SQ(4999.0) - 4.0); /* sqrt() is double, so store in temp FLOAT / x1 = (4999.0 - x)/2.0; x2 = (4999.0 + x)/2.0; (void) printf("x1 = %.16e x2 = %.16e\n",x1,x2); / roots computed directly */
if (x1 > x2) { xmax = x1; xmin = 1.0/x1; } else { xmax = x2; xmin = 1.0/x2; } (void) printf("xmin = %.16e xmax = %.16e\n",xmin,xmax); /* smaller from larger */
return 0; }
Output:
single-precision version x1 = 2.4414062500000000e-04 x2 = 4.9990000000000000e+ xmin = 2.0004001271445304e-04 xmax = 4.9990000000000000e+
double-precision version x1 = 2.0004001589768450e-04 x2 = 4.9989997999599836e+ xmin = 2.0004001600640290e-04 xmax = 4.9989997999599836e+
/* ** gmean.c */
#include #include
#define FLOAT float
#define NMAX 20
int main(void) { const FLOAT phi = (sqrt(5.0) - 1.0)/2.0; /* The Golden Mean */
** spin.c */
#include #include #include
#define N_MAX 2000 /* maximum array size */
#define SQ(x) ((x)*(x))
typedef struct { /* new type for storing particle data */ double r[3]; double v[3]; } DATA;
void get_inertia(DATA d,int n,float inertia[3][3]) { / computes inertia tensor of supplied data */
float r2; int i,j,k;
/* initialize */
for (j=0;j<3;j++) for (k=0;k<3;k++) inertia[j][k] = 0;
/* compute matrix */
for (i=0;i int main(int argc,char *argv[]) { extern float **convert_matrix(float *,long,long,long,long); extern void free_convert_matrix(float **,long,long,long,long); extern void gaussj(float **,int,float **,int);
DATA data[N_MAX]; double com_pos[3],com_vel[3],ang_mom[3]; float inertia[3][3],spin,tmp[3][1],a,b; int n;
if (argc != 2) { (void) fprintf(stderr,"Usage: %s dat-file\n",argv[0]); return 1; }
(void) printf("Reading...\n"); read_data(argv[1],data,&n,N_MAX); (void) printf("%i particle(s) read\n",n); assert(n > 0);
(void) printf("Analyzing...\n"); analyze_data(data,n,com_pos,com_vel,ang_mom); (void) printf("com pos = (%g,%g,%g)\ncom vel = (%g,%g,%g)\n" "ang mom = (%g,%g,%g)\n", com_pos[0],com_pos[1],com_pos[2], com_vel[0],com_vel[1],com_vel[2], ang_mom[0],ang_mom[1],ang_mom[2]);
(void) printf("Computing inertia tensor...\n"); get_inertia(data,n,inertia); (void) printf("inertia tensor:\n%g %g %g\n%g %g %g\n%g %g %g\n", inertia[0][0],inertia[0][1],inertia[0][2], inertia[1][0],inertia[1][1],inertia[1][2], inertia[2][0],inertia[2][1],inertia[2][2]);
(void) printf("Inverting...\n"); a = convert_matrix(&inertia[0][0],1,3,1,3); /* NRiC style matrix... / / gaussj expects its third argument to be a matrix... / tmp[0][0] = ang_mom[0]; tmp[1][0] = ang_mom[1]; tmp[2][0] = ang_mom[2]; b = convert_matrix(&tmp[0][0],1,3,1,1); gaussj(a,3,b,1); / the angular velocity is returned in b */
/* spin period is 2 pi over angular speed */
spin = 2*M_PI/sqrt(SQ(b[1][1]) + SQ(b[2][1]) + SQ(b[3][1]));
(void) printf("Spin = %g = %g h\n",spin,spin/3600.0);
/* free resources */
free_convert_matrix(b,1,1,1,3); free_convert_matrix(a,1,3,1,3);
return 0; }
Output:
Reading... 1908 particle(s) read Analyzing... com pos = (6.88422,-5.88037,-40.087) com vel = (0.000191741,-0.000238645,-0.000533925) ang mom = (-3077.62,-232.353,-923896) Computing inertia tensor... inertia tensor: 1.05501e+09 7.37927e+08 -1.31551e+ 7.37927e+08 1.73659e+09 - -1.31551e+07 -845270 2.27997e+ Inverting... Spin = 15496.5 = 4.30458 h
spin.sm
macro plot 3 { dev postencap spin$1.eps limits 0 1 0 1 relocate 0.5 1. putlabel 5 View along positive $1 axis limits $vmin $vmax $vmin $vmax box xlabel $2 (m) ylabel $3 (m) points $2 $ hardcopy }
data ps2.dat read {x 1 y 2 z 3}
