Commit 4d5ca304 authored by Reinhard Prix's avatar Reinhard Prix
Browse files

lalapps/pulsar: moved ephemeris-creation codes into new subdir CreateEphemeris/

   - Matt has ok'ed this
   - also removed all trailing whitespaces
   - passes 'make distcheck'
   - refs #20
Original: bf8375cf7a5b940c053df8fd0f91ab57936e1611
parent 8573808c
...@@ -33,6 +33,7 @@ AC_CONFIG_FILES([\ ...@@ -33,6 +33,7 @@ AC_CONFIG_FILES([\
src/pulsar/templateBanks/Makefile \ src/pulsar/templateBanks/Makefile \
src/pulsar/TimingTests/Makefile \ src/pulsar/TimingTests/Makefile \
src/pulsar/TDS_isolated/Makefile \ src/pulsar/TDS_isolated/Makefile \
src/pulsar/CreateEphemeris/Makefile \
src/inspiral/Makefile \ src/inspiral/Makefile \
src/inspiral/bayestar/Makefile \ src/inspiral/bayestar/Makefile \
src/inspiral/posterior/Makefile \ src/inspiral/posterior/Makefile \
......
...@@ -14,8 +14,8 @@ ...@@ -14,8 +14,8 @@
%% GNU General Public License for more details. %% GNU General Public License for more details.
%% %%
%% You should have received a copy of the GNU General Public License %% You should have received a copy of the GNU General Public License
%% along with with program; see the file COPYING. If not, write to the %% along with with program; see the file COPYING. If not, write to the
%% Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, %% Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
%% MA 02111-1307 USA %% MA 02111-1307 USA
%% %%
...@@ -26,12 +26,12 @@ ...@@ -26,12 +26,12 @@
%% which simplifies the expression quite a bit... %% which simplifies the expression quite a bit...
%% some useful constants %% some useful constants
YRSID_SI = 31558149.8; %% Sidereal year (1994), s YRSID_SI = 31558149.8; %% Sidereal year (1994), s
AU_SI = 1.4959787066e11; %% Astronomical unit, m AU_SI = 1.4959787066e11; %% Astronomical unit, m
C_SI = 299792458; %% Speed of light in vacuo, m s^-1 C_SI = 299792458; %% Speed of light in vacuo, m s^-1
## MLDC time-origin "t=0" (MUST be equal to setting in lal/packages/pulsar/include/LISAspecifics.h) ## MLDC time-origin "t=0" (MUST be equal to setting in lal/packages/pulsar/include/LISAspecifics.h)
LISA_TIME_ORIGIN = 700000000; LISA_TIME_ORIGIN = 700000000;
if ( nargin != 2 ) if ( nargin != 2 )
error ("\nNeed exactly two input arguments: start-GPS and end-GPS!"); error ("\nNeed exactly two input arguments: start-GPS and end-GPS!");
...@@ -42,7 +42,7 @@ endGPS = str2num ( argv{2} ); ...@@ -42,7 +42,7 @@ endGPS = str2num ( argv{2} );
duration = endGPS - startGPS; duration = endGPS - startGPS;
%% fields to write in the first line of ephem-file %% fields to write in the first line of ephem-file
gpsYr = ceil(startGPS); %% 'gpsYr' entry, not really used anywhere so we simply set it to start-time gpsYr = ceil(startGPS); %% 'gpsYr' entry, not really used anywhere so we simply set it to start-time
tStep = 14400; %% hardcoded default for now tStep = 14400; %% hardcoded default for now
nSteps = ceil( duration / tStep - 1e-6) + 1; nSteps = ceil( duration / tStep - 1e-6) + 1;
...@@ -52,7 +52,7 @@ lastGPS = startGPS + (nSteps-1) * tStep; ...@@ -52,7 +52,7 @@ lastGPS = startGPS + (nSteps-1) * tStep;
%% timesteps to compute LISA ephemeris for %% timesteps to compute LISA ephemeris for
tiGPS = [startGPS:tStep:lastGPS]' ; tiGPS = [startGPS:tStep:lastGPS]' ;
%% ----- implement LISA orbit Eq.(2.1) ----- %% ----- implement LISA orbit Eq.(2.1) -----
%% default at 't=0' %% default at 't=0'
kappa = 0; kappa = 0;
...@@ -83,7 +83,7 @@ if ( fid == -1 ) ...@@ -83,7 +83,7 @@ if ( fid == -1 )
error("Failed to open ephemeris-file '%s' for writing: error = %s", fname, msg ); error("Failed to open ephemeris-file '%s' for writing: error = %s", fname, msg );
endif endif
%% write header-line to file: %% write header-line to file:
fprintf (fid, "%d %f %d\n", gpsYr, tStep, nSteps ); fprintf (fid, "%d %f %d\n", gpsYr, tStep, nSteps );
%% write the time-series of [tGPS, pos, vel, acc ]: %% write the time-series of [tGPS, pos, vel, acc ]:
...@@ -92,6 +92,3 @@ ephem = [ tiGPS, xG, yG, zG, x1G, y1G, z1G, x2G, y2G, z2G ]; ...@@ -92,6 +92,3 @@ ephem = [ tiGPS, xG, yG, zG, x1G, y1G, z1G, x2G, y2G, z2G ];
fprintf (fid, "%.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g\n", ephem' ); fprintf (fid, "%.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g %.16g\n", ephem' );
fclose ( fid ); fclose ( fid );
## Process this file with automake to produce Makefile.in
AM_CPPFLAGS = -I$(top_srcdir)/src -I$(top_srcdir)/src/lalapps -I$(top_builddir)/src/lalapps -DPKG_DATA_DIR='"$(pkgdatadir)/"'
LDADD = $(top_builddir)/src/lalapps/liblalapps.la
bin_PROGRAMS = lalapps_create_solar_system_ephemeris \
lalapps_create_time_correction_ephemeris
lalapps_create_solar_system_ephemeris_SOURCES = create_solar_system_ephemeris.c
lalapps_create_time_correction_ephemeris_SOURCES = create_time_correction_ephemeris.c create_time_correction_ephemeris.h
...@@ -455,7 +455,7 @@ A[0], A[1], A[2]); ...@@ -455,7 +455,7 @@ A[0], A[1], A[2]);
gps = gps_2000 + (gps_JD[0] - jd_2000)*day + (INT4)(gps_JD[1]*day + 1.e-4); gps = gps_2000 + (gps_JD[0] - jd_2000)*day + (INT4)(gps_JD[1]*day + 1.e-4);
/* the int call and 1.d-4 are just to make sure gps is the ``right'' integer*/ /* the int call and 1.d-4 are just to make sure gps is the ``right'' integer*/
convert(gps_JD, time); convert(gps_JD, time);
pleph(coeffArray, time, inputs.target, R, fp); pleph(coeffArray, time, inputs.target, R, fp);
Rnow[0] = R[0]; Rnow[0] = R[0];
Rnow[1] = R[1]; Rnow[1] = R[1];
...@@ -463,7 +463,7 @@ A[0], A[1], A[2]); ...@@ -463,7 +463,7 @@ A[0], A[1], A[2]);
Vnow[0] = R[3]; Vnow[0] = R[3];
Vnow[1] = R[4]; Vnow[1] = R[4];
Vnow[2] = R[5]; Vnow[2] = R[5];
gps_JD[1] += halfinterval_jd; gps_JD[1] += halfinterval_jd;
if(gps_JD[1] >= 1.){ if(gps_JD[1] >= 1.){
gps_JD[1] -= 1.; gps_JD[1] -= 1.;
...@@ -576,7 +576,7 @@ INT4 fsizer(FILE *fp){ ...@@ -576,7 +576,7 @@ INT4 fsizer(FILE *fp){
/* flip bytes of ipt */ /* flip bytes of ipt */
endian_swap((CHAR*)&head1.data.ipt, sizeof(INT4), 36); endian_swap((CHAR*)&head1.data.ipt, sizeof(INT4), 36);
endian_swap((CHAR*)&head1.data.libratPtr, sizeof(INT4), 3); endian_swap((CHAR*)&head1.data.libratPtr, sizeof(INT4), 3);
/*** Calculate array size in the ephemeris */ /*** Calculate array size in the ephemeris */
for (i=0;i<12;i++){ for (i=0;i<12;i++){
if(head1.data.ipt[i][0] > kmx){ if(head1.data.ipt[i][0] > kmx){
...@@ -616,7 +616,7 @@ INT4 fsizer(FILE *fp){ ...@@ -616,7 +616,7 @@ INT4 fsizer(FILE *fp){
/* Initialise to ephemeris to the point at which the coefficient values start /* Initialise to ephemeris to the point at which the coefficient values start
*/ */
fseek(fp, 2*size*sizeof(REAL8), SEEK_SET); fseek(fp, 2*size*sizeof(REAL8), SEEK_SET);
/* flip bytes of values */ /* flip bytes of values */
endian_swap((CHAR*)&head1.data.au, sizeof(REAL8), 1); endian_swap((CHAR*)&head1.data.au, sizeof(REAL8), 1);
endian_swap((CHAR*)&head1.data.emrat, sizeof(REAL8), 1); endian_swap((CHAR*)&head1.data.emrat, sizeof(REAL8), 1);
......
...@@ -17,7 +17,7 @@ ...@@ -17,7 +17,7 @@
* MA 02111-1307 USA * MA 02111-1307 USA
*/ */
/* Many of these functions are (with occasional modification) taken directly /* Many of these functions are (with occasional modification) taken directly
* from the TEMPO2 software package http://www.atnf.csiro.au/research/pulsar/tempo2/ * from the TEMPO2 software package http://www.atnf.csiro.au/research/pulsar/tempo2/
* written by George Hobbs and Russell Edwards */ * written by George Hobbs and Russell Edwards */
...@@ -31,56 +31,56 @@ int verbose = 0; ...@@ -31,56 +31,56 @@ int verbose = 0;
int main(int argc, char **argv){ int main(int argc, char **argv){
inputParams_t inputParams; inputParams_t inputParams;
FILE *fp = NULL; FILE *fp = NULL;
int i=0, ne=0; int i=0, ne=0;
long double mjdtt=0.; long double mjdtt=0.;
double deltaT=0.; double deltaT=0.;
struct tm beg, fin; struct tm beg, fin;
char outputfile[MAXFNAME]; char outputfile[MAXFNAME];
/* get input arguments */ /* get input arguments */
get_input_args(&inputParams, argc, argv); get_input_args(&inputParams, argc, argv);
/* create output file name */ /* create output file name */
beg = *XLALGPSToUTC( &beg, (int)inputParams.startT ); beg = *XLALGPSToUTC( &beg, (int)inputParams.startT );
fin = *XLALGPSToUTC( &fin, (int)inputParams.endT ); fin = *XLALGPSToUTC( &fin, (int)inputParams.endT );
/* make sure the year in the filename only covers full years */ /* make sure the year in the filename only covers full years */
if ( beg.tm_yday != 0 && beg.tm_hour != 0 && beg.tm_min != 0 && beg.tm_sec != 0 ) beg.tm_year++; if ( beg.tm_yday != 0 && beg.tm_hour != 0 && beg.tm_min != 0 && beg.tm_sec != 0 ) beg.tm_year++;
if ( fin.tm_yday != 0 && fin.tm_hour != 0 && fin.tm_min != 0 && fin.tm_sec != 0 ) fin.tm_year--; if ( fin.tm_yday != 0 && fin.tm_hour != 0 && fin.tm_min != 0 && fin.tm_sec != 0 ) fin.tm_year--;
if( inputParams.et == TT2TCB ){ if( inputParams.et == TT2TCB ){
if( !sprintf(outputfile, "%s/te405_%d-%d.dat", if( !sprintf(outputfile, "%s/te405_%d-%d.dat",
inputParams.outputpath, beg.tm_year+1900, fin.tm_year+1900) ){ inputParams.outputpath, beg.tm_year+1900, fin.tm_year+1900) ){
fprintf(stderr, "Error... problem creating output file name!\n"); fprintf(stderr, "Error... problem creating output file name!\n");
exit(1); exit(1);
} }
} }
else if( inputParams.et == TT2TDB ){ else if( inputParams.et == TT2TDB ){
if( !sprintf(outputfile, "%s/tdb_%d-%d.dat", if( !sprintf(outputfile, "%s/tdb_%d-%d.dat",
inputParams.outputpath, beg.tm_year+1900, fin.tm_year+1900) ){ inputParams.outputpath, beg.tm_year+1900, fin.tm_year+1900) ){
fprintf(stderr, "Error... problem creating output file name!\n"); fprintf(stderr, "Error... problem creating output file name!\n");
exit(1); exit(1);
} }
} }
/* open the output file */ /* open the output file */
if( (fp = fopen(outputfile, "w")) == NULL ){ if( (fp = fopen(outputfile, "w")) == NULL ){
fprintf(stderr, "Error... can't open output file: %s!\n", outputfile); fprintf(stderr, "Error... can't open output file: %s!\n", outputfile);
exit(1); exit(1);
} }
ne = floor((inputParams.endT-inputParams.startT)/inputParams.interval); ne = floor((inputParams.endT-inputParams.startT)/inputParams.interval);
if( verbose ){ if( verbose ){
fprintf(stdout, "Creating ephemeris file from %lf to %lf, with %.2lf second\ fprintf(stdout, "Creating ephemeris file from %lf to %lf, with %.2lf second\
time steps\n", inputParams.startT, inputParams.endT, inputParams.interval); time steps\n", inputParams.startT, inputParams.endT, inputParams.interval);
} }
/* output header information on lines starting with a # comment */ /* output header information on lines starting with a # comment */
fprintf(fp, "# Build information for %s\n", argv[0]); fprintf(fp, "# Build information for %s\n", argv[0]);
fprintf(fp, "# Author: "LALAPPS_VCS_AUTHOR"\n"); fprintf(fp, "# Author: "LALAPPS_VCS_AUTHOR"\n");
...@@ -89,9 +89,9 @@ int main(int argc, char **argv){ ...@@ -89,9 +89,9 @@ int main(int argc, char **argv){
fprintf(fp, "#\n# Ephemeris creation command:-\n#\t"); fprintf(fp, "#\n# Ephemeris creation command:-\n#\t");
for( INT4 k=0; k<argc; k++ ) fprintf(fp, "%s ", argv[k]); for( INT4 k=0; k<argc; k++ ) fprintf(fp, "%s ", argv[k]);
fprintf(fp, "\n"); fprintf(fp, "\n");
CHAR *efile = strrchr(inputParams.ephemfile, '/'); CHAR *efile = strrchr(inputParams.ephemfile, '/');
fprintf(fp, "#\n# Time correction ephemeris file %s from TEMPO2 \ fprintf(fp, "#\n# Time correction ephemeris file %s from TEMPO2 \
(http://www.atnf.csiro.au/research/pulsar/tempo2/)\n", efile+1); (http://www.atnf.csiro.au/research/pulsar/tempo2/)\n", efile+1);
/* some information about the data */ /* some information about the data */
...@@ -103,17 +103,17 @@ int main(int argc, char **argv){ ...@@ -103,17 +103,17 @@ int main(int argc, char **argv){
fprintf(fp, "# with entries calculated at the Terrestrial Time in MJD via \ fprintf(fp, "# with entries calculated at the Terrestrial Time in MJD via \
the conversion:\n\ the conversion:\n\
# TT(MJD) = 44244 + (GPS + 51.184)/86400\n"); # TT(MJD) = 44244 + (GPS + 51.184)/86400\n");
fprintf(fp, "%lf\t%lf\t%lf\t%d\n", inputParams.startT, inputParams.endT, fprintf(fp, "%lf\t%lf\t%lf\t%d\n", inputParams.startT, inputParams.endT,
inputParams.interval, ne); inputParams.interval, ne);
for( i=0; i<ne; i++){ for( i=0; i<ne; i++){
/* convert GPS to equivalent TT value (in MJD) */ /* convert GPS to equivalent TT value (in MJD) */
mjdtt = MJDEPOCH + (inputParams.startT + (double)i*inputParams.interval + mjdtt = MJDEPOCH + (inputParams.startT + (double)i*inputParams.interval +
GPS2TT)/DAYSTOSEC; GPS2TT)/DAYSTOSEC;
/* these basically calculate the einstein delay */ /* these basically calculate the einstein delay */
/* using TEMPO2/TT to TCB/Teph file */ /* using TEMPO2/TT to TCB/Teph file */
if( inputParams.et == TT2TCB ) deltaT = IF_deltaT(mjdtt); if( inputParams.et == TT2TCB ) deltaT = IF_deltaT(mjdtt);
/* using TEMPO/TT to TDB file */ /* using TEMPO/TT to TDB file */
...@@ -124,16 +124,16 @@ the conversion:\n\ ...@@ -124,16 +124,16 @@ the conversion:\n\
fprintf(stderr, "Error... unknown ephemeris type!\n"); fprintf(stderr, "Error... unknown ephemeris type!\n");
exit(1); exit(1);
} }
/* print output */ /* print output */
fprintf(fp, "%.16lf\n", deltaT); fprintf(fp, "%.16lf\n", deltaT);
} }
IFTE_close_file(); IFTE_close_file();
fclose(fp); fclose(fp);
return 0; return 0;
} }
void get_input_args(inputParams_t *inputParams, int argc, char *argv[]){ void get_input_args(inputParams_t *inputParams, int argc, char *argv[]){
...@@ -153,7 +153,7 @@ void get_input_args(inputParams_t *inputParams, int argc, char *argv[]){ ...@@ -153,7 +153,7 @@ void get_input_args(inputParams_t *inputParams, int argc, char *argv[]){
char *program = argv[0]; char *program = argv[0];
char *tempo2path; char *tempo2path;
if(argc == 1){ if(argc == 1){
fprintf(stderr, "Error... no input parameters given! Try again!!!\n"); fprintf(stderr, "Error... no input parameters given! Try again!!!\n");
exit(0); exit(0);
...@@ -163,7 +163,7 @@ void get_input_args(inputParams_t *inputParams, int argc, char *argv[]){ ...@@ -163,7 +163,7 @@ void get_input_args(inputParams_t *inputParams, int argc, char *argv[]){
fprintf(stderr, "Error... TEMPO2 environment variable not set!\n"); fprintf(stderr, "Error... TEMPO2 environment variable not set!\n");
exit(1); exit(1);
} }
/* set defaults */ /* set defaults */
inputParams->interval = 60.; /* default to 60 seconds between */ inputParams->interval = 60.; /* default to 60 seconds between */
...@@ -171,7 +171,7 @@ void get_input_args(inputParams_t *inputParams, int argc, char *argv[]){ ...@@ -171,7 +171,7 @@ void get_input_args(inputParams_t *inputParams, int argc, char *argv[]){
inputParams->endT = 1072569615.; /* default end to 1 Jan 2014 00:00 UTC */ inputParams->endT = 1072569615.; /* default end to 1 Jan 2014 00:00 UTC */
inputParams->outputpath = NULL; inputParams->outputpath = NULL;
/* parse input arguments */ /* parse input arguments */
while ( 1 ){ while ( 1 ){
int option_index = 0; int option_index = 0;
...@@ -212,16 +212,16 @@ void get_input_args(inputParams_t *inputParams, int argc, char *argv[]){ ...@@ -212,16 +212,16 @@ void get_input_args(inputParams_t *inputParams, int argc, char *argv[]){
fprintf(stderr, "Unknown error while parsing options\n"); fprintf(stderr, "Unknown error while parsing options\n");
} }
} }
/* set the ephemeris file */ /* set the ephemeris file */
if( !strcmp(inputParams->ephemtype, "TEMPO") || if( !strcmp(inputParams->ephemtype, "TEMPO") ||
!strcmp(inputParams->ephemtype, "TDB") ){ /* using TEMPO/TDB file */ !strcmp(inputParams->ephemtype, "TDB") ){ /* using TEMPO/TDB file */
inputParams->et = TT2TDB; inputParams->et = TT2TDB;
sprintf(inputParams->ephemfile, "%s%s", tempo2path, TT2TDB_FILE); sprintf(inputParams->ephemfile, "%s%s", tempo2path, TT2TDB_FILE);
if( verbose ){ if( verbose ){
fprintf(stdout, "Using TEMPO-style TT to TDB conversion file: %s\n", fprintf(stdout, "Using TEMPO-style TT to TDB conversion file: %s\n",
inputParams->ephemfile); inputParams->ephemfile);
} }
} }
/* using TEMPO2/TCB/Teph file */ /* using TEMPO2/TCB/Teph file */
...@@ -229,12 +229,12 @@ void get_input_args(inputParams_t *inputParams, int argc, char *argv[]){ ...@@ -229,12 +229,12 @@ void get_input_args(inputParams_t *inputParams, int argc, char *argv[]){
!strcmp(inputParams->ephemtype, "TCB") ){ !strcmp(inputParams->ephemtype, "TCB") ){
inputParams->et = TT2TCB; inputParams->et = TT2TCB;
sprintf(inputParams->ephemfile, "%s%s", tempo2path, IFTEPH_FILE); sprintf(inputParams->ephemfile, "%s%s", tempo2path, IFTEPH_FILE);
if( verbose ){ if( verbose ){
fprintf(stdout, "Using TEMPO2-style TT to te405 conversion file: %s\n", fprintf(stdout, "Using TEMPO2-style TT to te405 conversion file: %s\n",
inputParams->ephemfile); inputParams->ephemfile);
} }
/* open and initialise file */ /* open and initialise file */
IFTE_init( inputParams->ephemfile ); IFTE_init( inputParams->ephemfile );
} }
...@@ -271,9 +271,9 @@ double FB_deltaT(long double mjd_tt, char fname[MAXFNAME]){ ...@@ -271,9 +271,9 @@ double FB_deltaT(long double mjd_tt, char fname[MAXFNAME]){
jda = read_double(); /* use jda as a dummy variable here */ jda = read_double(); /* use jda as a dummy variable here */
jda = read_double(); /* use jda as a dummy variable here */ jda = read_double(); /* use jda as a dummy variable here */
jda = read_double(); /* use jda as a dummy variable here */ jda = read_double(); /* use jda as a dummy variable here */
/* Use the corrected TT time and convert to Julian date */ /* Use the corrected TT time and convert to Julian date */
tdt = mjd_tt + 2400000.5; tdt = mjd_tt + 2400000.5;
if (tdt - (int)tdt >= 0.5) { if (tdt - (int)tdt >= 0.5) {
jda = (int)tdt + 0.5; jda = (int)tdt + 0.5;
jdb = tdt - (int)tdt - 0.5; jdb = tdt - (int)tdt - 0.5;
...@@ -282,7 +282,7 @@ double FB_deltaT(long double mjd_tt, char fname[MAXFNAME]){ ...@@ -282,7 +282,7 @@ double FB_deltaT(long double mjd_tt, char fname[MAXFNAME]){
jda = (int)tdt - 0.5; jda = (int)tdt - 0.5;
jdb = tdt - (int)tdt + 0.5; jdb = tdt - (int)tdt + 0.5;
} }
nr = (int)((jda-tdbd1)/tdbdt)+2; nr = (int)((jda-tdbd1)/tdbdt)+2;
if (nr < 1 || tdt > tdbd2){ if (nr < 1 || tdt > tdbd2){
printf("ERROR [CLK4]: Date %.10f out of range of TDB-TDT \ printf("ERROR [CLK4]: Date %.10f out of range of TDB-TDT \
table\n",(double)tdt); table\n",(double)tdt);
...@@ -297,30 +297,30 @@ table\n",(double)tdt); ...@@ -297,30 +297,30 @@ table\n",(double)tdt);
} }
t[0] = ((jda-((nr-2)*tdbdt+tdbd1))+jdb)/tdbdt; /* Fraction within record */ t[0] = ((jda-((nr-2)*tdbdt+tdbd1))+jdb)/tdbdt; /* Fraction within record */
t[1] = 1.0; /* Unused */ t[1] = 1.0; /* Unused */
/* Interpolation: call interp(buf,t,tdbncf,1, 1, 1, ctatv) */ /* Interpolation: call interp(buf,t,tdbncf,1, 1, 1, ctatv) */
np = 2; np = 2;
twot = 0.0; twot = 0.0;
pc[0] = 1.0; pc[1]=0.0; pc[0] = 1.0; pc[1]=0.0;
dna = 1.0; dna = 1.0;
dt1 = (int)(t[0]); dt1 = (int)(t[0]);
temp = dna * t[0]; temp = dna * t[0];
tc = 2.0*(fortran_mod(temp,1.0)+dt1)-1.0; tc = 2.0*(fortran_mod(temp,1.0)+dt1)-1.0;
if (tc != pc[1]){ if (tc != pc[1]){
np = 2; np = 2;
pc[1] = tc; pc[1] = tc;
twot = tc+tc; twot = tc+tc;
} }
if (np<tdbncf) { if (np<tdbncf) {
for (k=np+1;k<=tdbncf;k++) pc[k-1] = twot*pc[k-2]-pc[k-3]; for (k=np+1;k<=tdbncf;k++) pc[k-1] = twot*pc[k-2]-pc[k-3];
np = tdbncf; np = tdbncf;
} }
ctatv = 0.0; ctatv = 0.0;
for (j=tdbncf;j>=1;j--) ctatv = ctatv +pc[j-1]*buf[j-1]; for (j=tdbncf;j>=1;j--) ctatv = ctatv +pc[j-1]*buf[j-1];
close_file(); close_file();
return ctatv; return ctatv;
...@@ -339,17 +339,17 @@ int swapByte; ...@@ -339,17 +339,17 @@ int swapByte;
int open_file(char fname[MAXFNAME]){ int open_file(char fname[MAXFNAME]){
int i; int i;
/* NOTE MUST PUT BACK TO 0 FOR SOLARIS!!!!! */ /* NOTE MUST PUT BACK TO 0 FOR SOLARIS!!!!! */
union Convert{ union Convert{
char asChar[sizeof(int)]; char asChar[sizeof(int)];
int asInt; int asInt;
} intcnv; } intcnv;
intcnv.asInt = 1; intcnv.asInt = 1;
if (intcnv.asChar[0]==1) swapByte = 1; if (intcnv.asChar[0]==1) swapByte = 1;
else swapByte = 0; else swapByte = 0;
/* Look for blank character in filename */ /* Look for blank character in filename */
for (i=0;fname[i];i++) { for (i=0;fname[i];i++) {
if (fname[i]==' '){ if (fname[i]==' '){
...@@ -370,7 +370,7 @@ void close_file(void){ ...@@ -370,7 +370,7 @@ void close_file(void){
int read_int(void){ int read_int(void){
int i; int i;
union Convert{ union Convert{
char asChar[sizeof(int)]; char asChar[sizeof(int)];
int asInt; int asInt;
...@@ -421,7 +421,7 @@ void IFTE_init(const char fname[MAXFNAME]) { ...@@ -421,7 +421,7 @@ void IFTE_init(const char fname[MAXFNAME]) {
int ncon;