PROGRAM nbody C C This sample program writes an NBODY binary table extension using C the FITSIO library You should C be able to compile as follows: C C f77 -O -o nbody nbody.for -lfitsio C C (FITSIO: William Pence, HEASARC/NASA, V3.31 - but tested through C V3.415 - june 1994) C C The main program is designed as to hide most of the FITSIO specific C details, and merely show where the astronomy isdefined C C The input table is read from standard input. Within NEMO one C can make a sample table as follows: C C mkplummer tmp1 10 # create mass,pos,vel C hackforce tmp1 tmp2 # add phi and acc C snapdens tmp2 tmp3 # add aux C snapmask tmp3 tmp4 3:9 # add key C snapprint tmp4 m,x,y,z,vx,vy,vz,phi,ax,ay,az,aux,key > nbody.tab C C Conversion from this simple table to FITS is as follows: C C nbody < nbody.tab C C The output file, named 'nbody.fits', created by 'nbody' must not exist yet. C INTEGER MAXBODY PARAMETER (MAXBODY=1000) REAL mass(MAXBODY), pos(3,MAXBODY), vel(3,MAXBODY), * phi(MAXBODY), acc(3,MAXBODY) DOUBLE PRECISION aux(MAXBODY) INTEGER idx(MAXBODY), nbody, lun REAL time C -- get the data from a table nbody = MAXBODY CALL tab_get(nbody,mass,pos,vel,phi,acc,aux,idx) C -- open up file and initialize CALL nbf_open(lun,'nbody.fits') C -- write header and data seperately time = 0.0 CALL nbf_whead(lun, nbody, time) CALL nbf_wdata(lun, nbody, mass, pos, vel, phi, acc, aux, idx) C -- close up the shop CALL nbf_close(lun) END C*********************************************************************** SUBROUTINE tab_get(nbody,mass,pos,vel,phi,acc,aux,idx) INTEGER nbody REAL mass(*), pos(3,*), vel(3,*), phi(*), acc(3,*) DOUBLE PRECISION aux(*) INTEGER idx(*) C INTEGER i, maxbody C C -- read from stdin until no more C maxbody = nbody i = 1 100 CONTINUE READ(5,*,ERR=200,END=200) mass(i), * pos(1,i),pos(2,i),pos(3,i), * vel(1,i),vel(2,i),vel(3,i),phi(i), * acc(1,i),acc(2,i),acc(3,i),aux(i),idx(i) i=i+1 IF (i.LE.maxbody) GOTO 100 TYPE *,'Warning: bodies buffers exhausted' 200 CONTINUE nbody = i - 1 TYPE *,'Read ',nbody,' particles from input' END C*********************************************************************** SUBROUTINE nbf_open(lun,fname) INTEGER lun CHARACTER fname*(*) C- C Write the primary HDU - only needed once C- INTEGER status, bitpix, naxis, naxes(2), pcount, gcount LOGICAL simple, extend REAL version C C sanity check C CALL ftvers(version) TYPE *,'Using FITSIO Version ',version C C open the new FITS file C status=0 lun=21 CALL ftinit(lun,fname,2880,status) C C write a dummy primary header before the BINTABLE extension C this will only be a header, there will be no data C FITS requires this to be present. We're only wasting 2880 bytes C simple=.TRUE. bitpix=32 naxis=2 naxes(1)=0 naxes(2)=0 pcount=0 gcount=0 extend=.TRUE. C write the required primary array keywords CALL ftphpr(lun,simple,bitpix,naxis,naxes,pcount,gcount, * extend,status) C add today's date and some comments about the origin CALL ftpdat(lun,status) CALL ftpcom(lun, * 'This file was written by a sample program.', * status) CALL ftpkye(lun,'VERSION',version,3,'FITSIO version',status) CALL ftpcom(lun, * 'It was merely produced to show how to write an NBODY BINTABLE', * status) CALL ftpcom(lun,'Sample NEMO history: ',istatus) CALL ftphis(lun,'mkplummer t1 10',istatus) CALL ftphis(lun,'hackforce t1 t2',istatus) CALL ftphis(lun,'snapdens t2 t3',istatus) CALL ftphis(lun,'snapmask t3 t4 3:9',istatus) CALL ftphis(lun, * 'snapprint t4 m,x,y,z,vx,vy,vz,phi,ax,ay,az,aux,key > nbody.tab', * istatus) C define the primary array structure CALL ftpdef(lun,bitpix,naxis,naxes,pcount,gcount,status) C call error checker CALL nbf_error(status) END C*********************************************************************** SUBROUTINE nbf_whead(lun, nbody, time) INTEGER lun, nbody REAL time C- INTEGER MAXCOL PARAMETER (MAXCOL=16) INTEGER status INTEGER tfields,vardat,extver CHARACTER*8 tform(MAXCOL),ttype(MAXCOL),tunit(MAXCOL),extnam C C now create a binary table extension C C create a new empty extension CALL ftcrhd(lun,status) C define some global properties for this table C extnam='NBODY' extver=1 vardat=0 C define the columns (at this point you need to know how many, and C which columns to write. We've choosen not to to clutter the subroutine C argument list with this, but to make this routine general it C would be needed. C Note this in this case we didn't define tunit(), since UNITS='VIRIAL' C does it all. ttype(1)='MASS' tform(1)='1E' ttype(2)='POS' tform(2)='3E' ttype(3)='VEL' tform(3)='3E' ttype(4)='PHI' tform(4)='1E' ttype(5)='ACC' tform(5)='3E' ttype(6)='AUX' tform(6)='1D' ttype(7)='KEY' tform(7)='1J' tfields=7 C write the required keywords to the binary table extension CALL ftphbn(lun,nbody,tfields,ttype,tform,tunit,extnam, * vardat,status) CALL ftpdat(lun,status) CALL ftpkyj(lun,'EXTVER',extver, * 'NBODY Version (required)',status) CALL ftpkyj(lun,'NDIM',3, * 'Space dimensionality (required)',status) CALL ftpkye(lun,'TIME',time,9, * 'Snapshot time (optional [0])',status) CALL ftpkye(lun,'GRAVCON',1.0,1, * 'Gravitational constant (optional [1])',status) CALL ftpkys(lun,'UNITS','VIRIAL', * 'Unit Type (optional [VIRIAL])',status) CALL ftpkys(lun,'COORDSYS','CARTES', * 'Coordinate system (optional [CARTES])',status) C initialize the FITSIO parameters defining the structure of the table CALL ftbdef(lun,tfields,tform,vardat,nbody,status) C call error checker CALL nbf_error(status) END C*********************************************************************** SUBROUTINE nbf_wdata(lun,nbody,mass,pos,vel,phi,acc,aux,idx) INTEGER lun, nbody REAL mass(*), pos(3,*), vel(3,*), phi(*), acc(3,*) DOUBLE PRECISION aux(*) INTEGER idx(*) C- INTEGER status CALL ftpcle(lun,1,1,1,nbody,mass,status) CALL ftpcle(lun,2,1,1,nbody*3,pos,status) CALL ftpcle(lun,3,1,1,nbody*3,vel,status) CALL ftpcle(lun,4,1,1,nbody,phi,status) CALL ftpcle(lun,5,1,1,nbody*3,acc,status) CALL ftpcld(lun,6,1,1,nbody,aux,status) CALL ftpclj(lun,7,1,1,nbody,idx,status) C call error checker CALL nbf_error(status) END C*********************************************************************** SUBROUTINE nbf_close(lun) INTEGER lun C- INTEGER status C close the table CALL ftclos(lun,status) C call error checker CALL nbf_error(status) TYPE *,'Program appears to have run successfully' END C*********************************************************************** SUBROUTINE nbf_error(status) INTEGER status c CHARACTER errtxt*40 IF (status .GT. 0)THEN CALL ftgerr(status,errtxt) TYPE *,'*** ERROR - program did not run successfully ***' TYPE *,'status =',status,': ',errtxt STOP END IF END C***********************************************************************