Article 12757 of alt.sources: Path: nntpd.lkg.dec.com!pa.dec.com!crl.dec.com!crl.dec.com!bloom-beacon.mit.edu!spool.mu.edu!newspump.wustl.edu!crcnews.unl.edu!backbone!backbone!wayne From: wayne@backbone.uucp (Wayne Schlitt) Newsgroups: alt.sources Subject: XStar-2.0.0 Orbiting star system simulator Part05/06 Date: Sun, 27 Aug 1995 23:14:50 GMT Organization: The Backbone Cabal Lines: 2816 Sender: wayne@backbone (Wayne Schlitt) Message-ID: Reply-To: wayne@cse.unl.edu NNTP-Posting-Host: cse.unl.edu Originator: wayne@cse.unl.edu Submitted-by: wayne@backbone.UUCP Archive-name: XStar-2.0.0/part05 --- cut here --- # Continuation of Shell Archive, created by backbone!wayne # This is part 5 # To unpack the enclosed files, please use this file as input to the # Bourne (sh) shell. This can be most easily done by the command; # sh < thisfilename # ---------- file reset_hstep.c ---------- filename="reset_hstep.c" if [ -f $filename ] then echo File \"$filename\" already exists\! Skipping... filename=/dev/null # throw it away else echo extracting file reset_hstep.c... fi cat << 'END-OF-FILE' > $filename #include "xstar.h" #include "xstar_ext.h" /* ** XStar -- a animated n-body solver for X windows ** Copyright (C) 1995 Wayne Schlitt (wayne@cse.unl.edu) ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License as published by ** the Free Software Foundation; either version 2 of the License, or ** (at your option) any later version. */ /* * reset the hstep values to keep from overflowing the small hstep variable * in star_disp_type. */ void reset_hstep( star_disp_type *sd ) { int i; int min_hstep; int hstep_used; if( !num_disp_pt ) return; /* * find the minimum hstep */ min_hstep = sd->hsteps; for( i = 0; i < num_disp_pt; i++ ) { if( sd->disp_pts[i].star == DISP_PT_UNUSED ) continue; if( sd->disp_pts[i].hstep < min_hstep ) min_hstep = sd->disp_pts[i].hstep; } if( sd->erase_hstep < min_hstep ) min_hstep = sd->erase_hstep; /* make sure we won't have to call reset_hstep for at least 1k steps */ hstep_used = sd->hsteps - sd->erase_hstep; if( hstep_used > MAX_HSTEP - 1024 ) sd->erase_hstep += hstep_used - (MAX_HSTEP - 1024); /* * reset the hstep variables */ sd->hsteps -= min_hstep; sd->raw_hsteps -= min_hstep / sd->hstep_scale; sd->erase_hstep -= min_hstep; for( i = 0; i < num_disp_pt; i++ ) { if( sd->disp_pts[i].star == DISP_PT_UNUSED ) continue; sd->disp_pts[i].hstep -= min_hstep; } } END-OF-FILE if [ "$filename" != "/dev/null" ] then size=`wc -c < $filename` if [ $size != 1537 ] then echo $filename changed - should be 1537 bytes, not $size bytes fi chmod 660 $filename fi # ---------- file set_star_disp.c ---------- filename="set_star_disp.c" if [ -f $filename ] then echo File \"$filename\" already exists\! Skipping... filename=/dev/null # throw it away else echo extracting file set_star_disp.c... fi cat << 'END-OF-FILE' > $filename #include "xstar.h" #include "xstar_ext.h" /* ** XStar -- a animated n-body solver for X windows ** Copyright (C) 1995 Wayne Schlitt (wayne@cse.unl.edu) ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License as published by ** the Free Software Foundation; either version 2 of the License, or ** (at your option) any later version. */ /* * Set up the star display variables */ void set_star_disp( star_disp_type *sd, sys_param_type *sp ) { /* set up default erase_hstep adjust points */ sd->num_disp_pt_32 = num_disp_pt / 32; sd->num_disp_pt_16 = num_disp_pt / 16; sd->num_disp_pt_8 = num_disp_pt / 8; sd->num_disp_pt_4 = num_disp_pt / 4; /* handle special cases */ if( sp->do_bounce && sp->star_line ) sd->hsteps = sd->num_steps = MAX_HSTEP/2; else if( sp->do_bounce ) { if( sp->no_speed ) if( sp->num_collapsar == 0 ) sd->hsteps = sd->num_steps = 400 / sp->live_stars + 20; else sd->hsteps = sd->num_steps = 20; else sd->hsteps = sd->num_steps = 40 / sp->live_stars + 15; sd->num_disp_pt_32 = num_disp_pt * .875 + num_disp_pt / (32*4); sd->num_disp_pt_16 = num_disp_pt * .875 + num_disp_pt / (16*4); sd->num_disp_pt_8 = num_disp_pt * .875 + num_disp_pt / (8*4); sd->num_disp_pt_4 = num_disp_pt * .875 + num_disp_pt / (4*4); } else if( sp->few_stars ) { sd->hsteps = sd->num_steps = 15; sd->num_disp_pt_32 = num_disp_pt * .75 + num_disp_pt / (32*4); sd->num_disp_pt_16 = num_disp_pt * .75 + num_disp_pt / (16*4); sd->num_disp_pt_8 = num_disp_pt * .75 + num_disp_pt / (8*4); sd->num_disp_pt_4 = num_disp_pt * .75 + num_disp_pt / (4*4); } else if( sp->star_circle ) sd->hsteps = sd->num_steps = 400; else sd->hsteps = sd->num_steps = 200; /* if we don't have many display points, then let the tails grow slowly */ if( sd->num_disp_pt_4 - sd->num_disp_pt_8 < 1024 ) sd->num_disp_pt_4 += 1024 - (sd->num_disp_pt_4 - sd->num_disp_pt_8); /* set up other misc star_disp variables */ sd->erase_disps = 0; sd->live_erase = sp->live_stars; sd->next_free = 1; sd->next_disp = 1; sd->next_erase = 0; sd->erase_hstep = 0; sd->updt_hstep = sd->hsteps; sd->hstep_scale = fv_inv * (1./4); sd->raw_num_steps = sd->num_steps / sd->hstep_scale; sd->raw_hsteps = sd->hsteps / sd->hstep_scale; } void set_buffer_factor( star_disp_type *sd, sys_param_type *sp ) { int total_stars = sp->live_stars + sp->num_collapsar; sd->buffer_factor = max_disp_skip - total_stars*total_stars; if( sd->buffer_factor < 1 ) sd->buffer_factor = 1; } END-OF-FILE if [ "$filename" != "/dev/null" ] then size=`wc -c < $filename` if [ $size != 2692 ] then echo $filename changed - should be 2692 bytes, not $size bytes fi chmod 660 $filename fi # ---------- file move_stars_taylor3.c ---------- filename="move_stars_taylor3.c" if [ -f $filename ] then echo File \"$filename\" already exists\! Skipping... filename=/dev/null # throw it away else echo extracting file move_stars_taylor3.c... fi cat << 'END-OF-FILE' > $filename #include "xstar.h" #include "xstar_ext.h" #include /* ** XStar -- a animated n-body solver for X windows ** Copyright (C) 1995 Wayne Schlitt (wayne@cse.unl.edu) ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License as published by ** the Free Software Foundation; either version 2 of the License, or ** (at your option) any later version. */ /* * taylor series expansion to the 2nd order with the 3rd order approximated * by using the derivative the the 2nd degree Lagrange polynomial. */ static int num_rk4; void restart_taylor3( int n, double *m, point_2d *x, point_2d *v ) { num_rk4 = 3; } void init_taylor3( int n, double *m, point_2d *x, point_2d *v ) { init_rk4( n, m, x, v ); restart_taylor3( n, m, x, v ); } void move_stars_taylor3( point_2d *prev, point_2d *cur, double *m, point_2d *vk[], point_2d *ak[], int max_stars, sys_param_type *sp, star_disp_type *sd ) { double m_i, m_j; double previ_x, previ_y; double ai_x, ai_y; double Dt_a; /* derivative of the acceleration */ register int i, j, k; /* star index */ int tx, ty; int xpt, ypt; point_2d *v_0 = vk[0]; point_2d *v_1 = vk[1]; point_2d *a_0 = ak[0]; point_2d *a_1 = ak[1]; point_2d *a_2 = ak[2]; /* * we must use some other, self starting method for the first 3 * calls to the taylor method */ if( num_rk4 ) { --num_rk4; move_stars_rk4( prev, cur, m, vk, ak, max_stars, sp, sd ); return; } for( i = 0; i < max_stars; i++ ) { m_i = m[i]; if( m_i <= 0.0 ) continue; previ_x = prev[i].x; previ_y = prev[i].y; ai_x = a_0[i].x; ai_y = a_0[i].y; for( j = i+1; j < max_stars; j++ ) { double dx, dy, n2, f; register double f_d; m_j = m[j]; if( m_j <= 0.0 ) continue; dx = prev[j].x - previ_x; dy = prev[j].y - previ_y; n2 = dx*dx + dy*dy; if( n2 < collide_dist ) { a_0[i].x = ai_x; a_0[i].y = ai_y; collide( i, j, cur, prev, m, vk, ak, sp, sd ); m_i = m[i]; m_j = m[j]; previ_x = prev[i].x; previ_y = prev[i].y; ai_x = a_0[i].x; ai_y = a_0[i].y; } else { f = G / (sqrt(n2) * n2); f_d = dx * f; ai_x += f_d * m_j; a_0[j].x -= f_d * m_i; f_d = dy * f; ai_y += f_d * m_j; a_0[j].y -= f_d * m_i; } } /* * Move the star * * You want each step to be as small as possible and this is done * by making FV as large as possible. If the steps are small * enough, then we are approximating infinitesimal movements. * * The movement is based off of the formula: (with t=1) * x(t) = x0 + v t + 1/2 a t^2 + 1/(2*3) (Dt a) t^3 * v(t) = v0 + a t + 1/2 (Dt a) t^2 * * In theory, these formulas should include an infinite number * of derivatives of x(t), but the additional terms are not * easy to calculate. Even Dt a is just an estimate. * * * note: * I am currently using the derivative of the 2nd degree * Lagrange polynomial to estimate Dt_a. I tried using the * 3rd degree Lagrange Polynomial, but it didn't work as well. * Polynomials of high degree tend to not be very smooth so * their derivatives can make things worse. */ if( m_i >= collapsar ) { /* I think it looks better if the collapsars don't move ;-) */ v_1[i].x = v_1[i].y = 0.0; a_0[i].x = a_0[i].y = 0.0; ai_x = ai_y = 0.0; a_1[i].x = a_1[i].y = 0.0; xpt = cur[i].x = prev[i].x; ypt = cur[i].y = prev[i].y; } else { a_0[i].x = ai_x; /* Dt_a = (11./12)*ai_x - 1.5*a_1[i].x + .75*a_2[i].x - (1./6)*a_3[i].x; */ Dt_a = (3./4)*ai_x - a_1[i].x + (1./4)*a_2[i].x; v_0[i].x = v_1[i].x + ai_x + Dt_a; xpt = cur[i].x = previ_x + (v_1[i].x + ((1./2)*ai_x + (1./3)*Dt_a)); a_0[i].y = ai_y; /* Dt_a = (11./12)*ai_y - 1.5*a_1[i].y + .75*a_2[i].y - (1./6)*a_3[i].y; */ Dt_a = (3./4)*ai_y - a_1[i].y + (1./4)*a_2[i].y; v_0[i].y = v_1[i].y + ai_y + Dt_a; ypt = cur[i].y = previ_y + (v_1[i].y + ((1./2)*ai_y + (1./3)*Dt_a)); } ak[K-1][i].x = ak[K-1][i].y = 0; # include "plot_star.h" } } END-OF-FILE if [ "$filename" != "/dev/null" ] then size=`wc -c < $filename` if [ $size != 4277 ] then echo $filename changed - should be 4277 bytes, not $size bytes fi chmod 660 $filename fi # ---------- file move_stars_rk4.c ---------- filename="move_stars_rk4.c" if [ -f $filename ] then echo File \"$filename\" already exists\! Skipping... filename=/dev/null # throw it away else echo extracting file move_stars_rk4.c... fi cat << 'END-OF-FILE' > $filename #include "xstar.h" #include "xstar_ext.h" #include /* ** XStar -- a animated n-body solver for X windows ** Copyright (C) 1995 Wayne Schlitt (wayne@cse.unl.edu) ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License as published by ** the Free Software Foundation; either version 2 of the License, or ** (at your option) any later version. */ /* * Runge-Kutta method of 4th order */ static point_2d *ta_1, *ta_2, *ta_3; static int first_time = 1; void restart_rk4( int n, double *m, point_2d *x, point_2d *v ) { return; } void init_rk4( int n, double *m, point_2d *x, point_2d *v ) { int i; static int last_n = -1; default_init( n, m, x, v ); if( n != last_n ) { last_n = n; /* allocate storage */ if( !first_time ) { free( ta_1 ); free( ta_2 ); free( ta_3 ); } ta_1 = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); ta_2 = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); ta_3 = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); for( i = 0; i < n; i++ ) { ta_1[i].x = ta_1[i].y = 0; ta_2[i].x = ta_2[i].y = 0; ta_3[i].x = ta_3[i].y = 0; } } first_time = 0; } #define collide_rk4( i, j, cur, prev, m, vk, ak, ta_1, ta_2, ta_3, sp, sd ) do_collide_rk4( i, j, cur, prev, m, vk, ak, ta_1, ta_2, ta_3, sp, sd, __FILE__, __LINE__ ) static void do_collide_rk4( int i, int j, point_2d *cur, point_2d *prev, double *m, point_2d *vk[], point_2d *ak[], point_2d *ta_1, point_2d *ta_2, point_2d *ta_3, sys_param_type *sp, star_disp_type *sd, char *file, int line ) { /* collision: totally inelastic -> combine */ int k; double tmass; double orig_m_i = m[i]; num_collide++; if( verbose > 3 ) { prt_sys_const( max_stars, prev, m, vk[1] ); prt_sys_const_delta( orig_sc, max_stars, prev, m, vk[1] ); } tmass = 1. / (m[i] + m[j]); for( k = 0; k < K; k++ ) { vk[k][i].x = (m[i]*vk[k][i].x + m[j]*vk[k][j].x)*tmass; vk[k][i].y = (m[i]*vk[k][i].y + m[j]*vk[k][j].y)*tmass; } /* don't create a new collapsar */ if( m[i] >= collapsar ) m[i] += m[j]; else if ( m[j] >= collapsar ) { m[i] += m[j]; cur[i] = cur[j]; prev[i] = prev[j]; } else { /* the new location of the combined stars */ prev[i].x = (m[i]*prev[i].x + m[j]*prev[j].x)*tmass; prev[i].y = (m[i]*prev[i].y + m[j]*prev[j].y)*tmass; for( k = 0; k < K; k++ ) { ak[k][i].x = (m[i]*ak[k][i].x + m[j]*ak[k][j].x)*tmass; ak[k][i].y = (m[i]*ak[k][i].y + m[j]*ak[k][j].y)*tmass; } ta_1[i].x = (m[i]*ta_1[i].x + m[j]*ta_1[j].x)*tmass; ta_1[i].y = (m[i]*ta_1[i].y + m[j]*ta_1[j].y)*tmass; ta_2[i].x = (m[i]*ta_2[i].x + m[j]*ta_2[j].x)*tmass; ta_2[i].y = (m[i]*ta_2[i].y + m[j]*ta_2[j].y)*tmass; ta_3[i].x = (m[i]*ta_3[i].x + m[j]*ta_3[j].x)*tmass; ta_3[i].y = (m[i]*ta_3[i].y + m[j]*ta_3[j].y)*tmass; m[i] += m[j]; if( m[i] >= collapsar ) { /* This _should_, by construction, never happen. */ if( verbose > 0 ) fprintf( stderr, "%s:%d combined mass too large: m[%d]=%g\n", file, line, i, m[i]/fm ); m[i] = collapsar * .999; } } --sp->live_stars; sd->tstamp = time(NULL); set_buffer_factor( sd, sp ); if( verbose > 1 ) fprintf( stderr, "%s:%d live_stars=%d stars: %d,%d m: %g+%g=%g\n", file, line, sp->live_stars, i, j, orig_m_i/fm, m[j]/fm, m[i]/fm ); m[j] = 0.0; cur[j].x = cur[j].y = prev[j].x = prev[j].y = -1.0; for( k = 0; k < K; k++ ) { vk[k][j].x = vk[k][j].y = 0; ak[k][j].x = ak[k][j].y = 0; } ta_1[j].x = ta_1[j].y = 0; ta_2[j].x = ta_2[j].y = 0; ta_3[j].x = ta_3[j].y = 0; init_fna( max_stars, m, prev, vk[1] ); } void move_stars_rk4( point_2d *prev, point_2d *cur, double *m, point_2d *vk[], point_2d *ak[], int max_stars, sys_param_type *sp, star_disp_type *sd ) { double m_i, m_j; double previ_x, previ_y; double vi_x, vi_y; double a0i_x, a0i_y; double a1i_x, a1i_y; double ai_x, ai_y; double a_sum1, a_sum2; register int i, j, k; /* star index */ int tx, ty; int xpt, ypt; point_2d *v_0 = vk[0]; point_2d *v_1 = vk[1]; point_2d *a_0 = ak[0]; /* * initialize the ta_[1-3] arrays */ for( i = 0; i < max_stars; i++ ) { ta_1[i].x = ta_1[i].y = 0; ta_2[i].x = ta_2[i].y = 0; ta_3[i].x = ta_3[i].y = 0; } /* * move the stars */ for( i = 0; i < max_stars; i++ ) { m_i = m[i]; if( m_i <= 0.0 ) continue; previ_x = prev[i].x; previ_y = prev[i].y; /* * calculate the a_0 values */ ai_x = a_0[i].x; ai_y = a_0[i].y; for( j = i+1; j < max_stars; j++ ) { double dx, dy, n2, f; register double f_d; m_j = m[j]; if( m_j <= 0.0 ) continue; dx = prev[j].x - previ_x; dy = prev[j].y - previ_y; n2 = dx*dx + dy*dy; if( n2 < collide_dist ) { a_0[i].x = ai_x; a_0[i].y = ai_y; collide_rk4( i, j, cur, prev, m, vk, ak, ta_1, ta_2, ta_3, sp, sd ); m_i = m[i]; m_j = m[j]; previ_x = prev[i].x; previ_y = prev[i].y; ai_x = a_0[i].x; ai_y = a_0[i].y; } else { f = G / (sqrt(n2) * n2); f_d = dx * f; ai_x += f_d * m_j; a_0[j].x -= f_d * m_i; f_d = dy * f; ai_y += f_d * m_j; a_0[j].y -= f_d * m_i; } } a_0[i].x = ai_x; a_0[i].y = ai_y; /* * calculate the ta_1 values */ ai_x = ta_1[i].x; ai_y = ta_1[i].y; vi_x = v_1[i].x; vi_y = v_1[i].y; for( j = i+1; j < max_stars; j++ ) { double dx, dy, n2, f; register double f_d; m_j = m[j]; if( m_j <= 0.0 ) continue; dx = prev[j].x - previ_x + .5*(v_1[j].x - vi_x); dy = prev[j].y - previ_y + .5*(v_1[j].y - vi_y); n2 = dx*dx + dy*dy; if( n2 < collide_dist ) { ta_1[i].x = ai_x; ta_1[i].y = ai_y; collide_rk4( i, j, cur, prev, m, vk, ak, ta_1, ta_2, ta_3, sp, sd ); m_i = m[i]; m_j = m[j]; previ_x = prev[i].x; previ_y = prev[i].y; vi_x = v_1[i].x; vi_y = v_1[i].y; ai_x = ta_1[i].x; ai_y = ta_1[i].y; } else { f = G / (sqrt(n2) * n2); f_d = dx * f; ai_x += f_d * m_j; ta_1[j].x -= f_d * m_i; f_d = dy * f; ai_y += f_d * m_j; ta_1[j].y -= f_d * m_i; } } ta_1[i].x = ai_x; ta_1[i].y = ai_y; } for( i = 0; i < max_stars; i++ ) { m_i = m[i]; if( m_i <= 0.0 ) continue; previ_x = prev[i].x; previ_y = prev[i].y; /* * calculate the ta_2 values */ vi_x = v_1[i].x; vi_y = v_1[i].y; a0i_x = a_0[i].x; a0i_y = a_0[i].y; ai_x = ta_2[i].x; ai_y = ta_2[i].y; for( j = i+1; j < max_stars; j++ ) { double dx, dy, n2, f; register double f_d; m_j = m[j]; if( m_j <= 0.0 ) continue; dx = prev[j].x - previ_x + .5*(v_1[j].x - vi_x) + .25*(a_0[j].x - a0i_x); dy = prev[j].y - previ_y + .5*(v_1[j].y - vi_y) + .25*(a_0[j].y - a0i_y); n2 = dx*dx + dy*dy; if( n2 < collide_dist ) { ta_2[i].x = ai_x; ta_2[i].y = ai_y; collide_rk4( i, j, cur, prev, m, vk, ak, ta_1, ta_2, ta_3, sp, sd ); m_i = m[i]; m_j = m[j]; previ_x = prev[i].x; previ_y = prev[i].y; vi_x = v_1[i].x; vi_y = v_1[i].y; a0i_x = a_0[i].x; a0i_y = a_0[i].y; ai_x = ta_2[i].x; ai_y = ta_2[i].y; } else { f = G / (sqrt(n2) * n2); f_d = dx * f; ai_x += f_d * m_j; ta_2[j].x -= f_d * m_i; f_d = dy * f; ai_y += f_d * m_j; ta_2[j].y -= f_d * m_i; } } ta_2[i].x = ai_x; ta_2[i].y = ai_y; /* * calculate the ta_3 values */ a1i_x = ta_1[i].x; a1i_y = ta_1[i].y; ai_x = ta_3[i].x; ai_y = ta_3[i].y; for( j = i+1; j < max_stars; j++ ) { double dx, dy, n2, f; register double f_d; m_j = m[j]; if( m_j <= 0.0 ) continue; dx = prev[j].x - previ_x + v_1[j].x - vi_x + .5*(ta_1[j].x - a1i_x); dy = prev[j].y - previ_y + v_1[j].y - vi_y + .5*(ta_1[j].y - a1i_y); n2 = dx*dx + dy*dy; if( n2 < collide_dist ) { ta_3[i].x = ai_x; ta_3[i].y = ai_y; collide_rk4( i, j, cur, prev, m, vk, ak, ta_1, ta_2, ta_3, sp, sd ); m_i = m[i]; m_j = m[j]; previ_x = prev[i].x; previ_y = prev[i].y; vi_x = v_1[i].x; vi_y = v_1[i].y; a1i_x = ta_1[i].x; a1i_y = ta_1[i].y; ai_x = ta_3[i].x; ai_y = ta_3[i].y; } else { f = G / (sqrt(n2) * n2); f_d = dx * f; ai_x += f_d * m_j; ta_3[j].x -= f_d * m_i; f_d = dy * f; ai_y += f_d * m_j; ta_3[j].y -= f_d * m_i; } } ta_3[i].x = ai_x; ta_3[i].y = ai_y; /* * Move the star */ if( m_i >= collapsar ) { /* I think it looks better if the collapsars don't move ;-) */ v_1[i].x = v_1[i].y = 0.0; a_0[i].x = a_0[i].y = 0.0; ai_x = ai_y = 0.0; ta_1[i].x = ta_1[i].y = 0.0; xpt = cur[i].x = prev[i].x; ypt = cur[i].y = prev[i].y; } else { a_sum1 = ta_1[i].x + ta_2[i].x; a_sum2 = (1./6)*(a_0[i].x + a_sum1); xpt = cur[i].x = previ_x + (v_1[i].x + a_sum2); v_0[i].x = v_1[i].x + a_sum2 + (1./6)*(a_sum1 + ai_x); a_sum1 = ta_1[i].y + ta_2[i].y; a_sum2 = (1./6)*(a_0[i].y + a_sum1); ypt = cur[i].y = previ_y + (v_1[i].y + a_sum2); v_0[i].y = v_1[i].y + a_sum2 + (1./6)*(a_sum1 + ai_y); } ak[K-1][i].x = ak[K-1][i].y = 0; # include "plot_star.h" } } END-OF-FILE if [ "$filename" != "/dev/null" ] then size=`wc -c < $filename` if [ $size != 9550 ] then echo $filename changed - should be 9550 bytes, not $size bytes fi chmod 660 $filename fi # ---------- file move_stars_gpemce8.c ---------- filename="move_stars_gpemce8.c" if [ -f $filename ] then echo File \"$filename\" already exists\! Skipping... filename=/dev/null # throw it away else echo extracting file move_stars_gpemce8.c... fi cat << 'END-OF-FILE' > $filename #include "xstar.h" #include "xstar_ext.h" /* ** XStar -- a animated n-body solver for X windows ** Copyright (C) 1995 Wayne Schlitt (wayne@cse.unl.edu) ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License as published by ** the Free Software Foundation; either version 2 of the License, or ** (at your option) any later version. */ /* * Gragg's method of Polynomial-Extrapolation with Modifications for * Conservation of Energy. (i.e., this is a discrete mechanics method.) * * This method is from book _Numerical_Solutions_of_the_N-body_Problem_ by * Andrzej Marciniak. This appears to be his favorite method. */ #define H (1.) /* step size */ #define P (4) /* half the order of the approx */ #define MAX_SKIP (32) /* max # of correction skips */ /* * Notes: * * H should not be modified. To change the step size, modify fv instead. * * P is one half the order of the approximation. With P=4, FV=.25 you get * a very accurate approximation. With P=2, FV=.25, you get an approximation * that is slightly worse than the taylor approximation with FV=2, but the * taylor method is much faster. According to Andrzej's book, going beyond * P=4 won't get you much additional precision due to rounding errors. */ static int ibeta[P+1]; static int ibeta_2[P+1]; static double beta[P+1]; static double beta_2_inv[P+1]; static double beta_inv[P+1]; static double e0_2, a1, orig_a1, last_a1; static double alpha1[P+1], alpha2[P+1], bk[P+1], f[P], f1[P]; static point_2d *x1, *y, *z, *txk[P+1]; static point_2d *v1, *u, *w, *tvk[P+1]; static point_2d *sum2, *sum3; static double *sum, *sum1; static int first_time = 1; static int eps_skip; static int max_eps_skip; void restart_gpemce8( int n, double *m, point_2d *x, point_2d *v ) { return; } void init_gpemce8( int n, double *m, point_2d *x, point_2d *v ) { int i, j, k, l; double a, c; double sum4, sum5; static int last_n = -1; default_init( n, m, x, v ); if( n != last_n ) { last_n = n; /* allocate storage */ if( !first_time ) { free( x1 ); free( y ); free( z ); free( v1 ); free( u ); free( w ); free( sum ); free( sum1 ); free( sum2 ); free( sum3 ); for( k = 0; k < P+1; k++ ) { free( txk[k] ); free( tvk[k] ); } } x1 = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); y = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); z = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); v1 = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); u = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); w = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); sum = (double *)Safe_Malloc( sizeof( double ) * n ); sum1 = (double *)Safe_Malloc( sizeof( double ) * n ); sum2 = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); sum3 = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); for( k = 0; k < P+1; k++ ) { txk[k] = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); tvk[k] = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); } } /* calculate the total energy of the system (times 2) */ e0_2 = 0; for( i = 0; i < n; i++ ) { if( m[i] <= 0. ) continue; sum4 = 0; for( j = 0; j < n; j++ ) { if( m[j] <= 0. ) continue; if( i == j ) continue; sum4 += m[j] / sqrt( (x[i].x - x[j].x)*(x[i].x - x[j].x) + ( x[i].y - x[j].y )*( x[i].y - x[j].y )); } e0_2 += m[i]*( v[i].x*v[i].x + v[i].y*v[i].y - G*sum4 ); } if( first_time ) { /* step 1 */ ibeta[0] = 1; ibeta[1] = 2; ibeta[2] = 3; for( k = 3; k < P+1; k++ ) ibeta[k] = 2*ibeta[k-2]; for( k = 0; k < P+1; k++ ) { ibeta_2[k] = 2*ibeta[k]; beta[k] = ibeta[k]; beta_inv[k] = 1. / beta[k]; beta_2_inv[k] = .5 * beta_inv[k]; } /* step 2 */ for( k = 0; k < P+1; k++) { bk[k] = 1; for( l = 0; l < P+1; l++ ) { if( k == l ) continue; bk[k] *= ibeta[l]*ibeta[l]; } } for( i = 0; i < P; i++ ) f[i] = f1[i] = 1; for( l = 0; l < P-1; l++ ) { for( i = 0; i < P-1; i++ ) { if( l == i ) continue; a = bk[l] - bk[i]; f[l] *= (bk[P-1] - bk[i]) / a; f1[l] *= (bk[P] - bk[i]) / a; } f[l] *= bk[P-1]/bk[l]; f1[l] *= bk[P]/bk[l]; } a = c = 0; for( l = 0; l < P; l++ ) { a += f1[l]; c += f[l]; } a = 1 - a; c = 1 - c; alpha1[P] = 0; alpha1[P-1] = 1/c; alpha2[P] = 1; alpha2[P-1] = -a/c; for( l = 0; l < P-1; l++ ) { alpha1[l] = -f[l]/c; alpha2[l] = a*f[l]/c - f1[l]; } sum4 = sum5 = 0; for( k = 0; k < P+1; k++ ) { a = pow( beta[k], 2*P ); sum4 += alpha1[k]/a; sum5 += alpha2[k]/a; } orig_a1 = -sum4/sum5; } a1 = orig_a1; last_a1 = a1; eps_skip = -1; max_eps_skip = 4; first_time = 0; } void move_stars_gpemce8( point_2d *x, point_2d *x_new, double *m, point_2d *vk[], point_2d *ak[], int max_stars, sys_param_type *sp, star_disp_type *sd ) { double eps = 1./(1024. * 1024. * 1024. * 1024. ); int ibetak_2; double betak_2_inv; double betak_inv; double betak; int p1; int i, j, k; double ap, b, c, df, ff, r, r1, r2; int eps_calc; double dx, dy; double r_inv, r3_inv; double m_i; int tx, ty; int xpt, ypt; point_2d *v_new = vk[0]; point_2d *v = vk[1]; /* step 3 */ for( k = 0; k < P+1; k++ ) { ibetak_2 = ibeta_2[k]; betak = beta[k]; betak_inv = beta_inv[k]; betak_2_inv = beta_2_inv[k]; for( i = 0; i < max_stars; i++ ) { if( m[i] <= 0. ) continue; /* step 4 */ z[i] = x[i]; u[i] = v[i]; /* step 5 */ y[i].x = z[i].x + H*betak_2_inv*u[i].x; y[i].y = z[i].y + H*betak_2_inv*u[i].y; sum2[i].x = sum2[i].y = 0; } /* step 5 (continued) */ for( i = 0; i < max_stars; i++ ) { if( m[i] <= 0. ) continue; for( j = i+1; j < max_stars; j++ ) { if( m[j] <= 0. ) continue; dx = z[i].x - z[j].x; dy = z[i].y - z[j].y; r3_inv = 1 / (dx*dx + dy*dy); r3_inv = sqrt( r3_inv ) * r3_inv; sum2[i].x += m[j] * dx * r3_inv; sum2[j].x -= m[i] * dx * r3_inv; sum2[i].y += m[j] * dy * r3_inv; sum2[j].y -= m[i] * dy * r3_inv; } w[i].x = u[i].x - H*G*betak_2_inv*sum2[i].x; w[i].y = u[i].y - H*G*betak_2_inv*sum2[i].y; } for( p1 = 2; p1 < ibetak_2; p1++ ) { for( i = 0; i < max_stars; i++ ) { if( m[i] <= 0. ) continue; x1[i] = y[i]; v1[i] = w[i]; y[i].x = z[i].x + H*betak_inv*v1[i].x; y[i].y = z[i].y + H*betak_inv*v1[i].y; sum2[i].x = sum2[i].y = 0; } for( i = 0; i < max_stars; i++ ) { if( m[i] <= 0. ) continue; for( j = i+1; j < max_stars; j++ ) { if( m[j] <= 0. ) continue; dx = x1[i].x - x1[j].x; dy = x1[i].y - x1[j].y; r3_inv = 1 / (dx*dx + dy*dy); r3_inv = sqrt( r3_inv ) * r3_inv; sum2[i].x += m[j] * dx * r3_inv; sum2[j].x -= m[i] * dx * r3_inv; sum2[i].y += m[j] * dy * r3_inv; sum2[j].y -= m[i] * dy * r3_inv; } w[i].x = u[i].x - H*G*betak_inv*sum2[i].x; w[i].y = u[i].y - H*G*betak_inv*sum2[i].y; } for( i = 0; i < max_stars; i++ ) { z[i] = x1[i]; u[i] = v1[i]; } } /* step 6 */ for( i = 0; i < max_stars; i++ ) { sum2[i].x = sum2[i].y = 0; sum3[i].x = sum3[i].y = 0; } for( i = 0; i < max_stars; i++ ) { if( m[i] <= 0. ) continue; for( j = i+1; j < max_stars; j++ ) { if( m[j] <= 0. ) continue; dx = y[i].x - y[j].x; dy = y[i].y - y[j].y; r3_inv = 1 / (dx*dx + dy*dy); r3_inv = sqrt( r3_inv ) * r3_inv; sum2[i].x += m[j] * dx * r3_inv; sum2[j].x -= m[i] * dx * r3_inv; sum2[i].y += m[j] * dy * r3_inv; sum2[j].y -= m[i] * dy * r3_inv; dx = z[i].x - z[j].x + H*betak_inv*(w[i].x - w[j].x); dy = z[i].y - z[j].y + H*betak_inv*(w[i].y - w[j].y); r3_inv = 1 / (dx*dx + dy*dy); r3_inv = .5 * sqrt( r3_inv ) * r3_inv; sum3[i].x += m[j] * dx * r3_inv; sum3[j].x -= m[i] * dx * r3_inv; sum3[i].y += m[j] * dy * r3_inv; sum3[j].y -= m[i] * dy * r3_inv; } sum3[i].x += sum2[i].x; sum3[i].y += sum2[i].y; x1[i].x = (y[i].x + z[i].x)/2 + H*betak_2_inv*(w[i].x + u[i].x/2 - H*G*betak_2_inv*sum2[i].x ); x1[i].y = (y[i].y + z[i].y)/2 + H*betak_2_inv*(w[i].y + u[i].y/2 - H*G*betak_2_inv*sum2[i].y ); v1[i].x = (w[i].x + u[i].x)/2 - H*G*betak_2_inv*sum3[i].x; v1[i].y = (w[i].y + u[i].y)/2 - H*G*betak_2_inv*sum3[i].y; } for( i = 0; i < max_stars; i++ ) { if( m[i] <= 0. ) continue; txk[k][i] = x1[i]; tvk[k][i] = v1[i]; } } /* step 8 */ for( eps_calc = 0; eps_skip <= 1 && eps_calc < 16; eps_calc++ ) { ap = a1; ff = df = 0; for( i = 0; i < max_stars; i++ ) sum[i] = sum1[i] = 0; for( i = 0; i < max_stars; i++ ) { if( m[i] <= 0. ) continue; for( j = i+1; j < max_stars; j++ ) { if( m[j] <= 0. ) continue; r1 = r2 = 0; for( k = 0; k < P+1; k++ ) { b = txk[k][i].x - txk[k][j].x; r1 += (alpha1[k] + alpha2[k] * a1) * b; r2 += alpha2[k] * b; } r = r1*r1; c = r1 * r2; r1 = r2 = 0; for( k = 0; k < P+1; k++ ) { b = txk[k][i].y - txk[k][j].y; r1 += (alpha1[k] + alpha2[k] * a1) * b; r2 += alpha2[k] * b; } r += r1*r1; c += r1 * r2; r3_inv = 1 / r; r_inv = sqrt( r3_inv ); r3_inv *= c*r_inv; sum[i] += m[j] * r_inv; sum1[i] += m[j] * r3_inv; sum[j] += m[i] * r_inv; sum1[j] += m[i] * r3_inv; } r = c = 0; r1 = r2 = 0; for( k = 0; k < P+1; k++ ) { b = tvk[k][i].x; r1 += (alpha1[k] + alpha2[k] * a1) * b; r2 += alpha2[k] * b; } r += r1*r1; c += r1 * r2; r1 = r2 = 0; for( k = 0; k < P+1; k++ ) { b = tvk[k][i].y; r1 += (alpha1[k] + alpha2[k] * a1) * b; r2 += alpha2[k] * b; } r += r1*r1; c += r1 * r2; ff += m[i]*(r - G*sum[i]); df += m[i]*(2*c + G*sum1[i]); } a1 -= (ff - e0_2) / df; if( fabs( a1 - ap ) < eps ) break; } /* if( eps_calc > 0 || fabs( a1 - last_a1 ) > .125*eps ) */ if( eps_calc > 0 ) { last_a1 = a1; if( eps_skip >= 0 ) max_eps_skip = (max_eps_skip + 1) >> 1; eps_skip = -2; } /* * we don't need to update a1 _every_ time, so skip a few. updating a1 * twice in a row is important. It keeps the oscillations down. */ eps_skip++; if( eps_skip > max_eps_skip ) { eps_skip = 0; if( max_eps_skip < MAX_SKIP ) max_eps_skip += 1 + (max_eps_skip >> 3); } /* step 9 */ for( i = 0; i < max_stars; i++ ) { if( m[i] <= 0. || m[i] >= collapsar ) { x_new[i] = x[i]; v_new[i] = v[i]; } else { x_new[i].x = x_new[i].y = 0; v_new[i].x = v_new[i].y = 0; } } for( k = 0; k < P+1; k++ ) { b = alpha1[k] + alpha2[k] * a1; for( i = 0; i < max_stars; i++ ) { if( m[i] <= 0. || m[i] >= collapsar ) continue; x_new[i].x += b * txk[k][i].x; x_new[i].y += b * txk[k][i].y; v_new[i].x += b * tvk[k][i].x; v_new[i].y += b * tvk[k][i].y; } } /* * plot the points */ for( i = 0; i < max_stars; i++ ) { m_i = m[i]; if( m_i <= 0.0 ) continue; /* check for collisions */ for( j = i+1; j < max_stars; j++ ) { if( m[j] <= 0.0 ) continue; dx = x_new[i].x - x_new[j].x; dy = x_new[i].y - x_new[j].y; if( dx*dx + dy*dy < collide_dist ) collide( i, j, x, x_new, m, vk, ak, sp, sd ); } xpt = x_new[i].x; ypt = x_new[i].y; # define cur x_new # define prev x # define v_0 v # include "plot_star.h" } } END-OF-FILE if [ "$filename" != "/dev/null" ] then size=`wc -c < $filename` if [ $size != 12002 ] then echo $filename changed - should be 12002 bytes, not $size bytes fi chmod 660 $filename fi # ---------- file move_stars_ab7.c ---------- filename="move_stars_ab7.c" if [ -f $filename ] then echo File \"$filename\" already exists\! Skipping... filename=/dev/null # throw it away else echo extracting file move_stars_ab7.c... fi cat << 'END-OF-FILE' > $filename #include "xstar.h" #include "xstar_ext.h" #include /* ** XStar -- a animated n-body solver for X windows ** Copyright (C) 1995 Wayne Schlitt (wayne@cse.unl.edu) ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License as published by ** the Free Software Foundation; either version 2 of the License, or ** (at your option) any later version. */ /* * Adam-Bashford method of 7th order */ static int num_rk4; void restart_ab7( int n, double *m, point_2d *x, point_2d *v ) { num_rk4 = 7; } void init_ab7( int n, double *m, point_2d *x, point_2d *v ) { init_rk4( n, m, x, v ); restart_ab7( n, m, x, v ); } void move_stars_ab7( point_2d *prev, point_2d *cur, double *m, point_2d *vk[], point_2d *ak[], int max_stars, sys_param_type *sp, star_disp_type *sd ) { double m_i, m_j; double previ_x, previ_y; double ai_x, ai_y; register int i, j, k; /* star index */ int tx, ty; int xpt, ypt; point_2d *v_0 = vk[0]; point_2d *v_1 = vk[1]; point_2d *v_2 = vk[2]; point_2d *v_3 = vk[3]; point_2d *v_4 = vk[4]; point_2d *v_5 = vk[5]; point_2d *v_6 = vk[6]; point_2d *v_7 = vk[7]; point_2d *a_0 = ak[0]; point_2d *a_1 = ak[1]; point_2d *a_2 = ak[2]; point_2d *a_3 = ak[3]; point_2d *a_4 = ak[4]; point_2d *a_5 = ak[5]; point_2d *a_6 = ak[6]; /* * we must use some other, self starting method for the first 7 * calls to the Adam-Bashford 7th order method */ if( num_rk4 ) { --num_rk4; move_stars_rk4( prev, cur, m, vk, ak, max_stars, sp, sd ); return; } for( i = 0; i < max_stars; i++ ) { m_i = m[i]; if( m_i <= 0.0 ) continue; previ_x = prev[i].x; previ_y = prev[i].y; ai_x = a_0[i].x; ai_y = a_0[i].y; for( j = i+1; j < max_stars; j++ ) { double dx, dy, n2, f; register double f_d; m_j = m[j]; if( m_j <= 0.0 ) continue; dx = prev[j].x - previ_x; dy = prev[j].y - previ_y; n2 = dx*dx + dy*dy; if( n2 < collide_dist ) { a_0[i].x = ai_x; a_0[i].y = ai_y; collide( i, j, cur, prev, m, vk, ak, sp, sd ); m_i = m[i]; m_j = m[j]; previ_x = prev[i].x; previ_y = prev[i].y; ai_x = a_0[i].x; ai_y = a_0[i].y; } else { f = G / (sqrt(n2) * n2); f_d = dx * f; ai_x += f_d * m_j; a_0[j].x -= f_d * m_i; f_d = dy * f; ai_y += f_d * m_j; a_0[j].y -= f_d * m_i; } } /* * Move the star */ if( m_i >= collapsar ) { /* I think it looks better if the collapsars don't move ;-) */ v_1[i].x = v_1[i].y = 0.0; a_0[i].x = a_0[i].y = 0.0; ai_x = ai_y = 0.0; a_1[i].x = a_1[i].y = 0.0; xpt = cur[i].x = prev[i].x; ypt = cur[i].y = prev[i].y; } else { a_0[i].x = ai_x; v_0[i].x = v_1[i].x + (198721/60480.)*a_0[i].x - (447288/60480.)*a_1[i].x + (705549/60480.)*a_2[i].x - (688256/60480.)*a_3[i].x + (407139/60480.)*a_4[i].x - (134472/60480.)*a_5[i].x + (19087/60480.)*a_6[i].x; xpt = cur[i].x = prev[i].x + (198721/60480.)*v_1[i].x - (447288/60480.)*v_2[i].x + (705549/60480.)*v_3[i].x - (688256/60480.)*v_4[i].x + (407139/60480.)*v_5[i].x - (134472/60480.)*v_6[i].x + (19087/60480.)*v_7[i].x; a_0[i].y = ai_y; v_0[i].y = v_1[i].y + (198721/60480.)*a_0[i].y - (447288/60480.)*a_1[i].y + (705549/60480.)*a_2[i].y - (688256/60480.)*a_3[i].y + (407139/60480.)*a_4[i].y - (134472/60480.)*a_5[i].y + (19087/60480.)*a_6[i].y; ypt = cur[i].y = prev[i].y + (198721/60480.)*v_1[i].y - (447288/60480.)*v_2[i].y + (705549/60480.)*v_3[i].y - (688256/60480.)*v_4[i].y + (407139/60480.)*v_5[i].y - (134472/60480.)*v_6[i].y + (19087/60480.)*v_7[i].y; vk[K-1][i].x = vk[K-1][i].y = 0; } ak[K-1][i].x = ak[K-1][i].y = 0; # include "plot_star.h" } } END-OF-FILE if [ "$filename" != "/dev/null" ] then size=`wc -c < $filename` if [ $size != 3929 ] then echo $filename changed - should be 3929 bytes, not $size bytes fi chmod 660 $filename fi # ---------- file move_stars_am7.c ---------- filename="move_stars_am7.c" if [ -f $filename ] then echo File \"$filename\" already exists\! Skipping... filename=/dev/null # throw it away else echo extracting file move_stars_am7.c... fi cat << 'END-OF-FILE' > $filename #include "xstar.h" #include "xstar_ext.h" #include /* ** XStar -- a animated n-body solver for X windows ** Copyright (C) 1995 Wayne Schlitt (wayne@cse.unl.edu) ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License as published by ** the Free Software Foundation; either version 2 of the License, or ** (at your option) any later version. */ /* * 6th order Adam-Bashford predictor with a 7th order Adam-Moulton corrector, */ static point_2d *tloc, *tv, *ta; static int num_rk4; static int first_time = 1; void restart_am7( int n, double *m, point_2d *x, point_2d *v ) { num_rk4 = 5; } void init_am7( int n, double *m, point_2d *x, point_2d *v ) { int i; static int last_n = -1; init_rk4( n, m, x, v ); restart_am7( n, m, x, v ); if( n != last_n ) { last_n = n; /* allocate storage */ if( !first_time ) { free( tloc ); free( tv ); free( ta ); } tloc = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); tv = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); ta = (point_2d *)Safe_Malloc( sizeof( point_2d ) * n ); for( i = 0; i < n; i++ ) { tloc[i].x = tloc[i].y = 0; tv[i].x = tv[i].y = 0; ta[i].x = ta[i].y = 0; } } first_time = 0; } static void collide_am7( int i, int j, point_2d *cur, point_2d *prev, point_2d *tloc, double *m, point_2d *vk[], point_2d *ak[], point_2d *tv, point_2d *ta, sys_param_type *sp, star_disp_type *sd ) { /* collision: totally inelastic -> combine */ int k; double tmass; double orig_m_i = m[i]; num_collide++; if( verbose > 3 ) { prt_sys_const( max_stars, prev, m, vk[1] ); prt_sys_const_delta( orig_sc, max_stars, prev, m, vk[1] ); } tmass = 1. / (m[i] + m[j]); for( k = 0; k < K; k++ ) { vk[k][i].x = (m[i]*vk[k][i].x + m[j]*vk[k][j].x)*tmass; vk[k][i].y = (m[i]*vk[k][i].y + m[j]*vk[k][j].y)*tmass; } tv[i].x = (m[i]*tv[i].x + m[j]*tv[j].x)*tmass; tv[i].y = (m[i]*tv[i].y + m[j]*tv[j].y)*tmass; /* don't create a new collapsar */ if( m[i] >= collapsar ) m[i] += m[j]; else if ( m[j] >= collapsar ) { m[i] += m[j]; cur[i] = cur[j]; prev[i] = prev[j]; } else { /* the new location of the combined stars */ prev[i].x = (m[i]*prev[i].x + m[j]*prev[j].x)*tmass; prev[i].y = (m[i]*prev[i].y + m[j]*prev[j].y)*tmass; tloc[i].x = (m[i]*tloc[i].x + m[j]*tloc[j].x)*tmass; tloc[i].y = (m[i]*tloc[i].y + m[j]*tloc[j].y)*tmass; for( k = 0; k < K; k++ ) { ak[k][i].x = (m[i]*ak[k][i].x + m[j]*ak[k][j].x)*tmass; ak[k][i].y = (m[i]*ak[k][i].y + m[j]*ak[k][j].y)*tmass; } ta[i].x = (m[i]*ta[i].x + m[j]*ta[j].x)*tmass; ta[i].y = (m[i]*ta[i].y + m[j]*ta[j].y)*tmass; m[i] += m[j]; if( m[i] >= collapsar ) { /* This _should_, by construction, never happen. */ if( verbose > 0 ) fprintf( stderr, "%s:%d combined mass too large: m[%d]=%g\n", __FILE__, __LINE__, i, m[i]/fm ); m[i] = collapsar * .999; } } --sp->live_stars; sd->tstamp = time(NULL); set_buffer_factor( sd, sp ); if( verbose > 1 ) fprintf( stderr, "%s:%d live_stars=%d stars: %d,%d m: %g+%g=%g\n", __FILE__, __LINE__, sp->live_stars, i, j, orig_m_i/fm, m[j]/fm, m[i]/fm ); m[j] = 0.0; cur[j].x = cur[j].y = prev[j].x = prev[j].y = tloc[j].x = tloc[j].y = -1.0; for( k = 0; k < K; k++ ) { vk[k][j].x = vk[k][j].y = 0; ak[k][j].x = ak[k][j].y = 0; } tv[j].x = tv[j].y = 0; ta[j].x = ta[j].y = 0; init_fna( max_stars, m, tloc, tv ); } void move_stars_am7( point_2d *prev, point_2d *cur, double *m, point_2d *vk[], point_2d *ak[], int max_stars, sys_param_type *sp, star_disp_type *sd ) { double m_i, m_j; double previ_x, previ_y; double tloci_x, tloci_y; double ai_x, ai_y; register int i, j, k; /* star index */ int tx, ty; int xpt, ypt; point_2d *v_0 = vk[0]; point_2d *v_1 = vk[1]; point_2d *v_2 = vk[2]; point_2d *v_3 = vk[3]; point_2d *v_4 = vk[4]; point_2d *v_5 = vk[5]; point_2d *v_6 = vk[6]; point_2d *a_0 = ak[0]; point_2d *a_1 = ak[1]; point_2d *a_2 = ak[2]; point_2d *a_3 = ak[3]; point_2d *a_4 = ak[4]; point_2d *a_5 = ak[5]; /* * we must use some other, self starting method for the first 7 * calls to the Adam-Bashford predictor, Adam-Moulton corrector method */ if( num_rk4 ) { --num_rk4; move_stars_rk4( prev, cur, m, vk, ak, max_stars, sp, sd ); return; } /* * predictor step */ for( i = 0; i < max_stars; i++ ) { m_i = m[i]; if( m_i <= 0.0 ) continue; previ_x = prev[i].x; previ_y = prev[i].y; ai_x = a_0[i].x; ai_y = a_0[i].y; for( j = i+1; j < max_stars; j++ ) { double dx, dy, n2, f; register double f_d; m_j = m[j]; if( m_j <= 0.0 ) continue; dx = prev[j].x - previ_x; dy = prev[j].y - previ_y; n2 = dx*dx + dy*dy; if( n2 < collide_dist ) { a_0[i].x = ai_x; a_0[i].y = ai_y; collide( i, j, cur, prev, m, vk, ak, sp, sd ); m_i = m[i]; m_j = m[j]; previ_x = prev[i].x; previ_y = prev[i].y; ai_x = a_0[i].x; ai_y = a_0[i].y; } else { f = G / (sqrt(n2) * n2); f_d = dx * f; ai_x += f_d * m_j; a_0[j].x -= f_d * m_i; f_d = dy * f; ai_y += f_d * m_j; a_0[j].y -= f_d * m_i; } } /* * Move the star */ if( m_i >= collapsar ) { /* I think it looks better if the collapsars don't move ;-) */ v_1[i].x = v_1[i].y = 0.0; a_0[i].x = a_0[i].y = 0.0; ai_x = ai_y = 0.0; a_1[i].x = a_1[i].y = 0.0; tloc[i].x = prev[i].x; tloc[i].y = prev[i].y; } else { a_0[i].x = ai_x; a_0[i].y = ai_y; tv[i].x = v_1[i].x + (4227/1440.)*a_0[i].x - (7673/1440.)*a_1[i].x + (9482/1440.)*a_2[i].x - (6798/1440.)*a_3[i].x + (2627/1440.)*a_4[i].x - (425/1440.)*a_5[i].x; tv[i].y = v_1[i].y + (4227/1440.)*a_0[i].y - (7673/1440.)*a_1[i].y + (9482/1440.)*a_2[i].y - (6798/1440.)*a_3[i].y + (2627/1440.)*a_4[i].y - (425/1440.)*a_5[i].y; tloc[i].x = prev[i].x + (4227/1440.)*v_1[i].x - (7673/1440.)*v_2[i].x + (9482/1440.)*v_3[i].x - (6798/1440.)*v_4[i].x + (2627/1440.)*v_5[i].x - (425/1440.)*v_6[i].x; tloc[i].y = prev[i].y + (4227/1440.)*v_1[i].y - (7673/1440.)*v_2[i].y + (9482/1440.)*v_3[i].y - (6798/1440.)*v_4[i].y + (2627/1440.)*v_5[i].y - (425/1440.)*v_6[i].y; } } /* * corrector step */ for( i = 0; i < max_stars; i++ ) { m_i = m[i]; if( m_i <= 0.0 ) continue; tloci_x = tloc[i].x; tloci_y = tloc[i].y; ai_x = ta[i].x; ai_y = ta[i].y; for( j = i+1; j < max_stars; j++ ) { double dx, dy, n2, f; register double f_d; m_j = m[j]; if( m_j <= 0.0 ) continue; dx = tloc[j].x - tloci_x; dy = tloc[j].y - tloci_y; n2 = dx*dx + dy*dy; if( n2 < collide_dist ) { ta[i].x = ai_x; ta[i].y = ai_y; collide_am7( i, j, cur, prev, tloc, m, vk, ak, tv, ta, sp, sd ); m_i = m[i]; m_j = m[j]; tloci_x = tloc[i].x; tloci_y = tloc[i].y; ai_x = ta[i].x; ai_y = ta[i].y; } else { f = G / (sqrt(n2) * n2); f_d = dx * f; ai_x += f_d * m_j; ta[j].x -= f_d * m_i; f_d = dy * f; ai_y += f_d * m_j; ta[j].y -= f_d * m_i; } } /* * Move the star */ if( m_i >= collapsar ) { /* I think it looks better if the collapsars don't move ;-) */ v_1[i].x = v_1[i].y = 0.0; ta[i].x = ta[i].y = 0.0; ai_x = ai_y = 0.0; a_1[i].x = a_1[i].y = 0.0; xpt = cur[i].x = tloc[i].x; ypt = cur[i].y = tloc[i].y; } else { ta[i].x = ai_x; ta[i].y = ai_y; v_0[i].x = v_1[i].x + (19087/60480.)*ta[i].x + (65112/60480.)*a_0[i].x - (46461/60480.)*a_1[i].x + (37504/60480.)*a_2[i].x - (20211/60480.)*a_3[i].x + (6312/60480.)*a_4[i].x - (863/60480.)*a_5[i].x; v_0[i].y = v_1[i].y + (19087/60480.)*ta[i].y + (65112/60480.)*a_0[i].y - (46461/60480.)*a_1[i].y + (37504/60480.)*a_2[i].y - (20211/60480.)*a_3[i].y + (6312/60480.)*a_4[i].y - (863/60480.)*a_5[i].y; xpt = cur[i].x = prev[i].x + (19087/60480.)*tv[i].x + (65112/60480.)*v_1[i].x - (46461/60480.)*v_2[i].x + (37504/60480.)*v_3[i].x - (20211/60480.)*v_4[i].x + (6312/60480.)*v_5[i].x - (863/60480.)*v_6[i].x; ypt = cur[i].y = prev[i].y + (19087/60480.)*tv[i].y + (65112/60480.)*v_1[i].y - (46461/60480.)*v_2[i].y + (37504/60480.)*v_3[i].y - (20211/60480.)*v_4[i].y + (6312/60480.)*v_5[i].y - (863/60480.)*v_6[i].y; vk[K-1][i].x = vk[K-1][i].y = 0; } ta[i].x = ta[i].y = 0; ak[K-1][i].x = ak[K-1][i].y = 0; # include "plot_star.h" } } END-OF-FILE if [ "$filename" != "/dev/null" ] then size=`wc -c < $filename` if [ $size != 8867 ] then echo $filename changed - should be 8867 bytes, not $size bytes fi chmod 660 $filename fi # ---------- file dump_sys.c ---------- filename="dump_sys.c" if [ -f $filename ] then echo File \"$filename\" already exists\! Skipping... filename=/dev/null # throw it away else echo extracting file dump_sys.c... fi cat << 'END-OF-FILE' > $filename #include "xstar.h" #include "xstar_ext.h" /* ** XStar -- a animated n-body solver for X windows ** Copyright (C) 1995 Wayne Schlitt (wayne@cse.unl.edu) ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License as published by ** the Free Software Foundation; either version 2 of the License, or ** (at your option) any later version. */ /* * This routine dumps a star system into a format suitable for including * in another C program. This routine was used to create the init_sys_*() * routines. */ void dump_sys( point_2d *cur, point_2d *prev, double *m, point_2d *vk[], point_2d *ak[], int max_stars, sys_param_type *sp ) { int i, k; /* * dump the variables so that we can recreate this system later. */ fprintf( stderr, " sp->star_distrib = %d;\n", sp->star_distrib ); fprintf( stderr, " sp->size = %.16e;\n", sp->size ); fprintf( stderr, " sp->num_collapsar = %d;\n", sp->num_collapsar ); fprintf( stderr, " sp->mono_stars = %d;\n", sp->mono_stars ); fprintf( stderr, "\n" ); fprintf( stderr, " sp->star_circle = %d;\n", sp->star_circle ); fprintf( stderr, " sp->cir_dist = %d;\n", sp->cir_dist ); fprintf( stderr, "\n" ); fprintf( stderr, " sp->star_line = %d;\n", sp->star_line ); fprintf( stderr, " sp->rnd_spacing = %d;\n", sp->rnd_spacing ); fprintf( stderr, "\n" ); fprintf( stderr, " sp->do_bounce = %d;\n", sp->do_bounce ); fprintf( stderr, " sp->no_speed = %d;\n", sp->no_speed ); fprintf( stderr, " sp->few_stars = %d;\n", sp->few_stars ); fprintf( stderr, " sp->drift = %d;\n", sp->drift ); fprintf( stderr, "\n" ); fprintf( stderr, " sp->energy_fact = %d;\n", sp->energy_fact ); fprintf( stderr, " sp->calc_energy_fact = %d;\n", sp->calc_energy_fact ); fprintf( stderr, "\n" ); fprintf( stderr, " sp->min_stars = %d;\n", sp->min_stars ); fprintf( stderr, " sp->num_add = %d;\n", sp->num_add ); fprintf( stderr, " sp->live_stars = %d;\n", sp->live_stars ); for( i = 0; i < max_stars; i++ ) { if( m[i] <= 0 ) continue; fprintf( stderr, "\n/* star %d */\n", i ); fprintf( stderr, "\n m[%d] = %.16e;\n", i, m[i]/fm ); fprintf( stderr, " cur[%d].x = %.16e;\n", i, cur[i].x ); fprintf( stderr, " cur[%d].y = %.16e;\n", i, cur[i].y ); fprintf( stderr, " prev[%d].x = %.16e;\n", i, prev[i].x ); fprintf( stderr, " prev[%d].y = %.16e;\n", i, prev[i].y ); for( k = 0; k < K; k++ ) { if( vk[k][i].x != 0 || vk[k][i].y != 0 ) { fprintf( stderr, " vk[%d][%d].x = %.16e;\n", k, i, vk[k][i].x*fv ); fprintf( stderr, " vk[%d][%d].y = %.16e;\n", k, i, vk[k][i].y*fv ); } if( ak[k][i].x != 0 || ak[k][i].y != 0 ) { fprintf( stderr, " ak[%d][%d].x = %.16e;\n", k, i, ak[k][i].x*fv*fv ); fprintf( stderr, " ak[%d][%d].y = %.16e;\n", k, i, ak[k][i].y*fv*fv ); } } } } END-OF-FILE if [ "$filename" != "/dev/null" ] then size=`wc -c < $filename` if [ $size != 3186 ] then echo $filename changed - should be 3186 bytes, not $size bytes fi chmod 660 $filename fi # ---------- file init_sys_4.c ---------- filename="init_sys_4.c" if [ -f $filename ] then echo File \"$filename\" already exists\! Skipping... filename=/dev/null # throw it away else echo extracting file init_sys_4.c... fi cat << 'END-OF-FILE' > $filename #include "xstar.h" #include "xstar_ext.h" /* ** XStar -- a animated n-body solver for X windows ** Copyright (C) 1995 Wayne Schlitt (wayne@cse.unl.edu) ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License as published by ** the Free Software Foundation; either version 2 of the License, or ** (at your option) any later version. */ void init_sys_4( point_2d *cur, point_2d *prev, double *m, point_2d *vk[], point_2d *ak[], int max_stars, sys_param_type *sp ) { int i, k; sp->star_distrib = 0; sp->size = 2.1400000000000000e+02; sp->num_collapsar = 0; sp->mono_stars = 4; sp->star_circle = 0; sp->cir_dist = 0; sp->star_line = 0; sp->do_bounce = 0; sp->no_speed = 0; sp->few_stars = 0; sp->drift = 0; sp->min_stars = 1; sp->num_add = 0; sp->live_stars = 4; /* star 0 */ m[0] = 2.7865168539325840e+00; cur[0].x = 1.5835647927008775e+01; cur[0].y = 1.8923118463283608e+01; prev[0].x = 1.5835647927008775e+01; prev[0].y = 1.8923118463283608e+01; vk[0][0].x = 9.0941882306722933e-03; vk[0][0].y = 1.6257009059518187e-03; ak[0][0].x = 6.9931581889377990e-06; ak[0][0].y = -8.0395052303351502e-06; /* star 1 */ m[1] = 5.4831460674157304e+00; cur[1].x = 6.0521453125534116e+01; cur[1].y = 9.4644832984345086e+00; prev[1].x = 6.0521453125534116e+01; prev[1].y = 9.4644832984345086e+00; vk[0][1].x = 2.5826514208378845e-03; vk[0][1].y = -1.4348317268012971e-02; ak[0][1].x = -8.3617935918065052e-06; ak[0][1].y = -4.3824681501481231e-07; /* star 2 */ m[2] = 2.5168539325842696e+00; cur[2].x = 8.3340088264855723e+00; cur[2].y = -2.3248024535545774e+01; prev[2].x = 8.3340088264855723e+00; prev[2].y = -2.3248024535545774e+01; vk[0][2].x = -1.4385862343579154e-02; vk[0][2].y = -1.3381252605483343e-03; ak[0][2].x = 3.0889420909050185e-06; ak[0][2].y = 9.3058808330760419e-06; /* star 3 */ m[3] = 5.2134831460674160e+00; cur[3].x = -7.6139068509249114e+01; cur[3].y = -8.8449218374312917e+00; prev[3].x = -7.6139068509249114e+01; prev[3].y = -8.8449218374312917e+00; vk[0][3].x = -6.3202458968509645e-04; vk[0][3].y = 1.4867554009648933e-02; ak[0][3].x = 3.5653642534100792e-06; ak[0][3].y = 2.6539749189938063e-07; /* scale things to the current accuracy */ for( i = 0; i < max_stars; i++ ) { if( m[i] <= 0 ) continue; m[i] *= fm; for( k = 0; k < K; k++ ) { vk[k][i].x *= fv_inv; vk[k][i].y *= fv_inv; ak[k][i].x *= fv_inv*fv_inv; ak[k][i].y *= fv_inv*fv_inv; } } } END-OF-FILE if [ "$filename" != "/dev/null" ] then size=`wc -c < $filename` if [ $size != 2871 ] then echo $filename changed - should be 2871 bytes, not $size bytes fi chmod 660 $filename fi # ---------- file init_sys_8.c ---------- filename="init_sys_8.c" if [ -f $filename ] then echo File \"$filename\" already exists\! Skipping... filename=/dev/null # throw it away else echo extracting file init_sys_8.c... fi cat << 'END-OF-FILE' > $filename #include "xstar.h" #include "xstar_ext.h" /* ** XStar -- a animated n-body solver for X windows ** Copyright (C) 1995 Wayne Schlitt (wayne@cse.unl.edu) ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License as published by ** the Free Software Foundation; either version 2 of the License, or ** (at your option) any later version. */ void init_sys_8( point_2d *cur, point_2d *prev, double *m, point_2d *vk[], point_2d *ak[], int max_stars, sys_param_type *sp ) { int i, k; sp->star_distrib = 0; sp->size = 3.0300000000000000e+02; sp->num_collapsar = 0; sp->mono_stars = 8; sp->star_circle = 0; sp->cir_dist = 0; sp->star_line = 0; sp->do_bounce = 0; sp->no_speed = 0; sp->few_stars = 0; sp->drift = 0; sp->min_stars = 1; sp->num_add = 0; sp->live_stars = 8; /* star 0 */ m[0] = 1.1468208092485550e+00; cur[0].x = -2.7038883602489072e+00; cur[0].y = 6.4089177018335889e+01; prev[0].x = -2.7038883602489072e+00; prev[0].y = 6.4089177018335889e+01; vk[0][0].x = 1.3226751051071546e-02; vk[0][0].y = -2.2643948440514610e-03; ak[0][0].x = -2.1461637019132762e-06; ak[0][0].y = -3.8329229061867005e-06; /* star 1 */ m[1] = 2.2566473988439308e+00; cur[1].x = 6.0566200308691172e+01; cur[1].y = 5.0696810406423374e+01; prev[1].x = 6.0566200308691172e+01; prev[1].y = 5.0696810406423374e+01; vk[0][1].x = 1.1326416168617813e-02; vk[0][1].y = -1.2983888267850539e-02; ak[0][1].x = -3.7978097624060328e-06; ak[0][1].y = -4.9601977809216730e-06; /* star 2 */ m[2] = 1.0358381502890175e+00; cur[2].x = -1.3325368021270076e+01; cur[2].y = 4.3795680059746633e+00; prev[2].x = -1.3325368021270076e+01; prev[2].y = 4.3795680059746633e+00; vk[0][2].x = -1.4657477777562395e-03; vk[0][2].y = 8.3098025452551811e-03; ak[0][2].x = 1.4837373960419059e-06; ak[0][2].y = 3.2846997147130252e-06; /* star 3 */ m[3] = 2.7560693641618497e+00; cur[3].x = 5.0296086642545227e+01; cur[3].y = 3.3756236820588086e+00; prev[3].x = 5.0296086642545227e+01; prev[3].y = 3.3756236820588086e+00; vk[0][3].x = -3.7597216177549879e-03; vk[0][3].y = -8.3720019727285024e-03; ak[0][3].x = -1.9157531920650712e-06; ak[0][3].y = 3.8658856402276140e-06; /* star 4 */ m[4] = 2.1271676300578037e+00; cur[4].x = 6.8270243736874789e+01; cur[4].y = -1.1595121023145049e+02; prev[4].x = 6.8270243736874789e+01; prev[4].y = -1.1595121023145049e+02; vk[0][4].x = -1.4273624489641538e-02; vk[0][4].y = -8.7326800300004132e-03; ak[0][4].x = -9.7523640637911302e-07; ak[0][4].y = 1.7221688018067123e-06; /* star 5 */ m[5] = 2.0531791907514454e+00; cur[5].x = -6.1922907126650045e+01; cur[5].y = 8.2905813493962910e+01; prev[5].x = -6.1922907126650045e+01; prev[5].y = 8.2905813493962910e+01; vk[0][5].x = 1.4057051300353198e-02; vk[0][5].y = 9.7559665347273592e-03; ak[0][5].x = 4.3248908536138410e-06; ak[0][5].y = -6.5951618951828798e-06; /* star 6 */ m[6] = 1.9236994219653181e+00; cur[6].x = -4.8033093344416187e+01; cur[6].y = 4.7578502331255180e+01; prev[6].x = -4.8033093344416187e+01; prev[6].y = 4.7578502331255180e+01; vk[0][6].x = 3.8568338976656485e-03; vk[0][6].y = 9.4050873297924203e-03; ak[0][6].x = 2.2103324433326556e-06; ak[0][6].y = 3.6631219712243281e-06; /* star 7 */ m[7] = 2.7005780346820809e+00; cur[7].x = -6.8161132128329925e+01; cur[7].y = -8.0295095516064634e+01; prev[7].x = -6.8161132128329925e+01; prev[7].y = -8.0295095516064634e+01; vk[0][7].x = -1.2873818798120700e-02; vk[0][7].y = 9.9296098370418636e-03; ak[0][7].x = 1.3764904099135169e-06; ak[0][7].y = 1.6155758780025582e-06; /* scale things to the current accuracy */ for( i = 0; i < max_stars; i++ ) { if( m[i] <= 0 ) continue; m[i] *= fm; for( k = 0; k < K; k++ ) { vk[k][i].x *= fv_inv; vk[k][i].y *= fv_inv; ak[k][i].x *= fv_inv*fv_inv; ak[k][i].y *= fv_inv*fv_inv; } } } END-OF-FILE if [ "$filename" != "/dev/null" ] then size=`wc -c < $filename` if [ $size != 4459 ] then echo $filename changed - should be 4459 bytes, not $size bytes fi chmod 660 $filename fi # ---------- file init_sys_8b.c ---------- filename="init_sys_8b.c" if [ -f $filename ] then echo File \"$filename\" already exists\! Skipping... filename=/dev/null # throw it away else echo extracting file init_sys_8b.c... fi cat << 'END-OF-FILE' > $filename #include "xstar.h" #include "xstar_ext.h" /* ** XStar -- a animated n-body solver for X windows ** Copyright (C) 1995 Wayne Schlitt (wayne@cse.unl.edu) ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License as published by ** the Free Software Foundation; either version 2 of the License, or ** (at your option) any later version. */ void init_sys_8b( point_2d *cur, point_2d *prev, double *m, point_2d *vk[], point_2d *ak[], int max_stars, sys_param_type *sp ) { int i, k; sp->star_distrib = 0; sp->size = 3.0300000000000000e+02; sp->num_collapsar = 0; sp->mono_stars = 8; sp->star_circle = 0; sp->cir_dist = 0; sp->star_line = 0; sp->do_bounce = 0; sp->no_speed = 0; sp->few_stars = 0; sp->drift = 0; sp->min_stars = 1; sp->num_add = 0; sp->live_stars = 8; /* star 0 */ m[0] = 3.2820512820512824e+00; cur[0].x = -5.9426018201326272e+01; cur[0].y = 1.7647796986427664e+01; prev[0].x = -5.9426018201326272e+01; prev[0].y = 1.7647796986427664e+01; vk[0][0].x = 4.8740542740093326e-03; vk[0][0].y = 1.0900511132255125e-02; ak[0][0].x = 2.5173820214687024e-06; ak[0][0].y = -2.4452551992655210e-06; /* star 1 */ m[1] = 2.0968660968660968e+00; cur[1].x = 9.3599198133270978e-01; cur[1].y = -4.0668403188234521e+01; prev[1].x = 9.3599198133270978e-01; prev[1].y = -4.0668403188234521e+01; vk[0][1].x = -7.7175472328441120e-03; vk[0][1].y = -3.7383560620019456e-03; ak[0][1].x = -3.2319296217630353e-06; ak[0][1].y = 1.7991971056596420e-06; /* star 2 */ m[2] = 2.0512820512820515e+00; cur[2].x = -5.9384119374432700e+01; cur[2].y = 9.8236339712784101e+01; prev[2].x = -5.9384119374432700e+01; prev[2].y = 9.8236339712784101e+01; vk[0][2].x = 1.5656891177842101e-02; vk[0][2].y = 8.0608082916574265e-03; ak[0][2].x = 8.2457296929936326e-07; ak[0][2].y = -3.3782755206993198e-06; /* star 3 */ m[3] = 9.3447293447293456e-01; cur[3].x = -7.0405575139421671e+01; cur[3].y = -2.7001090804229296e+01; prev[3].x = -7.0405575139421671e+01; prev[3].y = -2.7001090804229296e+01; vk[0][3].x = -7.6525422837526603e-03; vk[0][3].y = 1.6937754974762298e-02; ak[0][3].x = 6.5073079181275713e-06; ak[0][3].y = 3.6029649981805933e-06; /* star 4 */ m[4] = 2.2336182336182335e+00; cur[4].x = 9.2302766801034323e+01; cur[4].y = -6.0862441826543247e+01; prev[4].x = 9.2302766801034323e+01; prev[4].y = -6.0862441826543247e+01; vk[0][4].x = -1.0083930782485945e-02; vk[0][4].y = -1.3905687652693736e-02; ak[0][4].x = -2.1482127913348859e-06; ak[0][4].y = 2.1385211861947785e-06; /* star 5 */ m[5] = 7.5213675213675213e-01; cur[5].x = -7.2473741939590468e+00; cur[5].y = 1.6015571178554193e+01; prev[5].x = -7.2473741939590468e+00; prev[5].y = 1.6015571178554193e+01; vk[0][5].x = 8.7707293778696745e-03; vk[0][5].y = -9.8131165679121637e-04; ak[0][5].x = -4.0382422553598351e-06; ak[0][5].y = -3.0253523635815457e-06; /* star 6 */ m[6] = 2.9173789173789175e+00; cur[6].x = 8.7255538788646888e+01; cur[6].y = 2.6925933985080562e+01; prev[6].x = 8.7255538788646888e+01; prev[6].y = 2.6925933985080562e+01; vk[0][6].x = 5.5994912063467676e-03; vk[0][6].y = -1.4944457539263970e-02; ak[0][6].x = -2.0118695310450290e-06; ak[0][6].y = -1.6850985387665771e-06; /* star 7 */ m[7] = 1.7321937321937322e+00; cur[7].x = -4.3062922104262000e+01; cur[7].y = -5.9796442419559682e+01; prev[7].x = -4.3062922104262000e+01; prev[7].y = -5.9796442419559682e+01; vk[0][7].x = -1.4541568021132618e-02; vk[0][7].y = 8.7153380803802905e-03; ak[0][7].x = 2.5674930287103935e-06; ak[0][7].y = 5.9061580476882484e-06; /* scale things to the current accuracy */ for( i = 0; i < max_stars; i++ ) { if( m[i] <= 0 ) continue; m[i] *= fm; for( k = 0; k < K; k++ ) { vk[k][i].x *= fv_inv; vk[k][i].y *= fv_inv; ak[k][i].x *= fv_inv*fv_inv; ak[k][i].y *= fv_inv*fv_inv; } } } END-OF-FILE if [ "$filename" != "/dev/null" ] then size=`wc -c < $filename` if [ $size != 4465 ] then echo $filename changed - should be 4465 bytes, not $size bytes fi chmod 660 $filename fi echo end of this archive file.... exit 0 --- cut here --- -- Wayne Schlitt can not assert the truth of all statements in this article and still be consistent.