!___________________________________________________________________________ ! HI SURVEYS SOFTWARE PACKAGE !___________________________________________________________________________ ! ! This file contains FORTRAN software for use in extracting spectra ! from the HI surveys of Heiles and Habing (HH), Parkes (PK), ! Argentina (AR), Bell Labs (BL), and Weaver and Williams (WW). ! ! By using these programs and the data that they access, the user ! agrees to the following: ! 1. If data from any survey are used in a publication, the user ! will reference the original authors of the survey, as if the ! original authors had themselves provided the data to the user, ! by quoting the appropriate publication in the usual fashion. ! 2. Whenever data obtained from the NSSDC are used in a publication ! acknowledgment shall be made to Dr. Wayne H. Warren Jr. and/or ! the Astronomical Data Center/National Space Science Data Center ! for storing and distributing the data. ! 3. If the survey data or software described herein are used, ! acknowledgment shall be made to William T. Reach for creating ! the tape and writing the software. ! ! There is no implied guarantee that these routines will perform the ! tasks they claim. In other words, user beware--check the results ! before using them. ! ! The programs included are, in order of appearance: ! DUMP_NSDC, which copies data from tape to disk, ! SPECTRUM_NSDC, which extracts individual spectra from the disk file, ! and the subroutine packages ! DECODE_NSDC, which contains routines to decode the spectra, and ! TAPE HANDLING ROUTINES. ! ! Note: The tape handling routines will not work on any machine other than ! a VAX under the VMS operating system; they are included only for ! completeness. ! program dump_nsdc ! copyright (c) 1989 by William T. Reach, University of California, Berkeley ! ! use this program to copy data from tape to disk ! the tape will be positioned to the beginning of the desired survey, ! based on the assumption that the order of datasets on the tape ! is HH, PK, AR, BL (for the "high-latitude tape") ! (if your tape is ordered differently, change the tape positioning section) ! character block*16384,record*4096 integer spectrum(512) integer nspectra(5),irecord_length(5),nchannels(5) data nspectra / 10, 10, 20, 20, 10/ data irecord_length / 643, 420, 396, 703,1204/ data nchannels / 100, 64, 56, 124, 238/ 100 write (*,102) 102 format (/,' NSSDC HI SURVEY TAPE DUMPER',/, . x,' (1) Heiles-Habing',/, . x,' (2) Parkes',/, . x,' (3) Argentina',/, . x,' (4) Bell Labs',/, . x,' (5) Weaver-Williams',/, . x,' Enter desired survey: ',$) read *,isurvey write (*,'(x,a,$)') 'Which tape drive: ' read *,idrive call getchannel (ichan,idrive) ! position tape at beginning of desired survey write (*,'(x,a,$)') 'Does the tape need to be positioned ("y"/n) ? ' read (*,'(a)') yn if (yn.eq.'n' .or. yn.eq.'N') goto 200 if (isurvey.ge.1 .and. isurvey.le.4) then print *,'Tape labelled "high-latitude" should be mounted' else if (isurvey.eq.5) then print *,'Tape labelled "Weaver-Williams" should be mounted' else print *,'invalid survey...try again' goto 100 endif call rewind_tape (ichan) if (isurvey.gt.1 .and. isurvey.lt.5) call skipfile(ichan,isurvey-1) ! position tape at first desired block ! (obtain information from the table of galactic latitude vs. block number) 200 print *,'first block containing desired data:' read *,iblock1 print *,'number of blocks to read (0=all):' read *,nblock_out if (iblock1.ne.1) call skipblock (ichan,iblock1-1) nblock = iblock1-1 write (*,250) 250 format (/,' writing to output file HI; 1 spectrum per record...',/) ! copy the data from tape to disk iok = 0 nspec = 0 ib = 0 do while (iok.eq.0) ib = ib + 1 if (nblock_out.ne.0 .and. ib.gt.nblock_out) goto 900 nbreq = 16384 nbytein = 0 call readblock (ichan,block,nbreq,nbytein) if (nbytein.eq.0) then print *,'end of file detected...' goto 900 endif nblock = nblock + 1 print *,'read block ',nblock,' # bytes = ',nbytein if (ib.eq.1) then irecl = irecord_length(isurvey) open (unit=1,file='hi',status='new',recordtype='fixed', . form='formatted',recl=irecl) endif ns = nbytein /irecl do i = 1,ns record(1:irecl) = block((i-1)*irecl+1:i*irecl) write (1,'(a)') record(1:irecl) nspec = nspec + 1 enddo if (mod(ib,10).eq.0) print *,'finished block ',ib enddo close (1) 900 print *,'done...wrote ',nspec,' spectra' end program spectrum_nsdc c copyright (c) 1989 by William T. Reach, University of California, Berkeley ! ! This routine searches through the HI survey file on disk. It finds the ! observation nearest a desired position, to within DIST_MAX degrees. ! The spectrum at this position is written to a three-column file ! containing the channel number, VLSR of the channel, and antenna ! temperature, respectively. ! ! The input file is assumed to be assigned to logical name HI; this file ! should have been produced by the DUMP_NSDC program. ! ! The output file is the default unit 1 file, FOR001 ! character record*4096 real spectrum(512),vlsr(512) integer nspectra(5),irecord_length(5),nchannels(5) data nspectra / 10, 10, 20, 20, 10/ data irecord_length / 643, 420, 396, 703,1204/ data nchannels / 100, 64, 56, 124, 238/ real dist_max(5) data dist_max /0.5,1.0,1.0,1.1,0.4/ 100 write (*,102) 102 format (/,' NSSDC HI SURVEY SPECTRUM FINDER',/, . x,' (1) Heiles-Habing',/, . x,' (2) Parkes',/, . x,' (3) Argentina',/, . x,' (4) Bell Labs',/, . x,' (5) Weaver-Williams',/, . x,' Enter desired survey: ',$) read *,isurvey irecl=irecord_length(isurvey) open (unit=20,file='hi',status='old',readonly, . form='formatted',recordtype='fixed',recl=irecl) print *,'desired position (gl,gb):' read *,gl0,gb0 ifound_it = 0 nrec = 0 do while (ifound_it.eq.0) nrec = nrec + 1 read (20,'(a)',end=800) record(1:irecl) call decode_spectrum (record,isurvey,gl,gb,ra,dec, . vlsr,spectrum) dist = sqrt((gl-gl0)**2 + (gb-gb0)**2) if (dist.lt.dist_max(isurvey)) then ifound_it = 1 print *,'I found it...' print *,'gl = ',gl,' gb = ',gb do i = 1,nchannels(isurvey) write (1,210) i,vlsr(i),spectrum(i) enddo endif enddo stop 210 format (x,i4,2(x,g)) 800 print *,'couldn''t find it...' end !___________________________________________________________________________ ! DECODING ROUTINES ! ! The following routines decode individual spectra of the various HI surveys, ! returning the antenna temperature as a function of VLSR. ! ! The first routine, DECODE_SPECTRUM, is a driver for the others and is the ! only one called by higher-level routines. ! ! INPUT VARIABLES: ! record: character buffer containing the record from the disk file ! isurvey: the survey number, 1=HH, 2=PK, 3=AR, 4=BL, 5=WW ! ! OUTPUT VARIABLES: ! gl,gb: galactic coordinates of the observed position (decimal degrees) ! ra,dec: equatorial coordinates, decimal hours, degrees ! vlsr: array containing the VLSR (km/s) of each spectral channel ! spectrum: array containing the antenna temperature of each spectral channel ! subroutine decode_spectrum (record,isurvey,gl,gb,ra,dec,vlsr,spectrum) ! copyright (c) 1989 by William T. Reach, University of California, Berkeley character record*4096 real spectrum(512),vlsr(512) integer nspectra(5),irecord_length(5),nchannels(5) data nspectra / 10, 10, 20, 20, 10/ data irecord_length / 643, 420, 396, 703,1204/ data nchannels / 100, 64, 56, 124, 238/ if (isurvey.eq.1) then call decode_hh (record,isurvey,gl,gb,ra,dec,vlsr,spectrum) else if (isurvey.eq.2) then call decode_pk (record,isurvey,gl,gb,ra,dec,vlsr,spectrum) else if (isurvey.eq.3) then call decode_ar (record,isurvey,gl,gb,ra,dec,vlsr,spectrum) else if (isurvey.eq.4) then call decode_bl (record,isurvey,gl,gb,ra,dec,vlsr,spectrum) else if (isurvey.eq.5) then call decode_ww (record,isurvey,gl,gb,ra,dec,vlsr,spectrum) endif end ! subroutine decode_hh (record,isurvey,gl,gb,ra,dec,vlsr,spectrum) ! copyright (c) 1989 by William T. Reach, University of California, Berkeley character record*4096 real spectrum(512),vlsr(512) integer*4 ispec(100),igl,igb read (record(1:643),901) igl,igb, . irah,iram,iras,idecd,idecm,idecs, . idatey,idatem,idated,ilst,icvel,istat, . ispec 901 format (2i6,3i2,i3,2i2,3i2,2i5,i1,1x,100i6) gl = real(igl)/1000. gb = real(igb)/1000. ra = real(irah) + real(iram)/60. + real(iras)/3600. dec = abs(real(idecd)) + real(idecm)/60. + real(idecs)/3600. if (idecd.lt.0) dec = -dec date = real(idatey) + real(idatem-1)/12. + real(idated-1)/(30.6*12.) cvel = real(icvel)/1000. do i = 1,100 spectrum(i) = real(ispec(i))/100. enddo vwide = (30.e3/1420.4058e6)*2.998e5 vnarr = (10.e3/1420.4058e6)*2.998e5 if (cvel.eq.0.) then do i = 1,8 vlsr(i) = real(i-1)*vwide - 90.76 enddo do i = 9,95 vlsr(i) = real(i-9)*vnarr/2. - 45.38 enddo do i = 96,100 vlsr(i) = real(i-96)*vwide + 46.43 enddo else if (cvel.eq.-8.44) then do i = 1,8 vlsr(i) = real(i-1)*vwide - 82.32 enddo do i = 9,95 vlsr(i) = real(i-9)*vnarr/2. - 36.94 enddo do i = 96,100 vlsr(i) = real(i-96)*vwide + 54.87 enddo else print *,'DECODE_HH> unrecognized central velocity' endif end ! subroutine decode_pk (record,isurvey,gl,gb,ra,dec,vlsr,spectrum) ! copyright (c) 1989 by William T. Reach, University of California, Berkeley character record*4096 real spectrum(512),vlsr(512) integer*4 ispec(64),igl,igb,icv read (record(1:420),901) igl,igb, . irah,iram,iras,idecd,idecm,idecs, . icv,istat, . ispec 901 format (2i6,3i2,i3,2i2,i7,i3,1x,64i6) gl = real(igl)/1000. gb = real(igb)/1000. ra = real(irah) + real(iram)/60. + real(iras)/3600. dec = abs(real(idecd)) + real(idecm)/60. + real(idecs)/3600. if (idecd.lt.0) dec = -dec cvel = real(icv)/1000. do i = 1,64 spectrum(i) = real(ispec(i))/100. vlsr(i) = cvel + real(i-1)*6.99 enddo end ! subroutine decode_ar (record,isurvey,gl,gb,ra,dec,vlsr,spectrum) ! copyright (c) 1989 by William T. Reach, University of California, Berkeley character record*4096 real spectrum(512),vlsr(512) integer*4 ispec(56),ibadscan(5),igl,igb,iv1,iv56 read (record(1:396),901) igl,igb, . irah,iram,iras,idecd,idecm,idecs, . ilsth,ilstm,ilsts,idatey,idatem,idated, . (ibadscan(k),k=1,5), . iv1,iv56,istat, . ispec 901 format (2i6,3i2,i3,2i2,6i2,5i2,2i6,i1,56i6) gl = real(igl)/1000. gb = real(igb)/1000. ra = real(irah) + real(iram)/60. + real(iras)/3600. dec = abs(real(idecd)) + real(idecm)/60. + real(idecs)/3600. if (idecd.lt.0) dec = -dec date = real(idatey) + real(idatem-1)/12. + real(idated-1)/(30.6*12.) v1 = real(iv1)/100. v56 = real(iv56)/100. dv = (v56-v1)/55. do i = 1,56 spectrum(i) = real(ispec(i))/100. vlsr(i) = v1 + (i-1)*dv enddo end ! subroutine decode_bl (record,isurvey,gl,gb,ra,dec,vlsr,spectrum) ! copyright (c) 1989 by William T. Reach, University of California, Berkeley character record*4096 real spectrum(512),vlsr(512) integer*4 ispec(124),igl,igb,icv integer*4 icoeffn(3),icoeffw(3) real coeffn(3),coeffw(3) read (record(1:703),901) igl,igb, . irah,iram,iras,idecd,idecm,idecs, . ich0,icv,itint,irms,nchbase, . icoeffn,icoeffw, . ispec 901 format (2i6,3i2,i3,2i2,i4,i6,i4,i5,i3,6i6,124i5) gl = real(igl)/100. gb = real(igb)/100. ra = real(irah) + real(iram)/60. + real(iras)/3600. dec = abs(real(idecd)) + real(idecm)/60. + real(idecs)/3600. if (idecd.lt.0) dec = -dec ch0 = real(ich0)/10. cvel = real(icv)/100. do i = 1,3 coeffn(i) = real(icoeffn(i))/1.e4 coeffw(i) = real(icoeffw(i))/1.e4 enddo do i = 1,124 spectrum(i) = real(ispec(i))/100. vlsr(i) = cvel + (real(i)-ch0)*5.2767 enddo end ! subroutine decode_ww (record,isurvey,gl,gb,ra,dec,vlsr,spectrum) ! copyright (c) 1989 by William T. Reach, University of California, Berkeley character record*4096 real spectrum(512),vlsr(512) integer*4 ispec(238),icv,igl,igb read (record(1:1204),901) igl,igb,icv,ispec 901 format (i4,i5,i4,i1,238i5) gl = real(igl)/10. gb = real(igb)/100. cvel = real(icv)/10. do i = 1,238 spectrum(i) = real(ispec(i))/10. vlsr(i) = cvel + (real(i)-122)*1.0553 enddo return do i = 1,238 spectrum(i) = 0. vlsr(i) = 0. enddo end !______________________________________________________________________________ ! TAPE HANDLING ROUTINES ! ! The following routines are a tape driver package that perform the interface ! between FORTRAN and the VAX/VMS system I/O routines. They will not work on ! any other operating system. They assume the tape drives are assigned to ! the names MUB0: and MUB1:, and that there are two tape drives. If your ! system is configured differently, change the routine GETCHANNEL accordingly. ! c*****************************************************************GETCHANNEL c c ASSIGN CHANNEL TO THE TAPE DRIVE c N.B. ALL I/O OPERATIONS TO TAPE DRIVE REQUIRE A CHANNEL c subroutine getchannel (channel,drive) ! copyright (c) 1989 by William T. Reach, University of California, Berkeley integer*4 status,sys$assign,channel,drive if (drive .eq. -1) then status = sys$assign ('MUB0:',channel,,) if (.not. status) status = sys$assign ('MUB1:',channel,,) if (.not. status) call lib$signal(%val(status)) return endif if (drive .gt. 2 .or. drive .lt. -1) then print *,'tape drive ',drive,' does not exist .. access denied' return endif if (drive .eq. 0) status = sys$assign ('MUB0:',channel,,) if (drive .eq. 1) status = sys$assign ('MUB1:',channel,,) if (.not. status) call lib$signal(%val(status)) end c c******************************************************************READBLOCK c c READ NEXT LOGICAL BLOCK ON TAPE c BUFFER = CHARACTER VARIABLE TO RECIEVE BLOCK c SIZE = NUMBER OF BYTES TO RECEIVE c NBYTES = NUMBER OF BYTES ACTUALLY READ INTO BUFFER c subroutine readblock (channel,buffer,size,nbytes) ! copyright (c) 1989 by William T. Reach, University of California, Berkeley integer*4 status,sys$qiow,iosb(2),channel,size integer*2 func character*(*) buffer func = 33 status = sys$qiow(, %val(channel), %val(func), iosb, , , + %ref(buffer(1:1)), %val(size), , , , ) if (.not. status) call lib$signal(%val(status)) do 10 n = 30,12,-1 if (iosb(1) .lt. 2**n) goto 10 iosb(1) = iosb(1) - 2**n nbytes = nbytes + 2**(n-16) 10 continue return end c c*****************************************************************REWIND_TAPE c c REWIND TAPE c subroutine rewind_tape (channel) ! copyright (c) 1989 by William T. Reach, University of California, Berkeley integer*4 status,sys$qiow,iosb(2),channel integer*2 func func = 36 status = sys$qiow(,%val(channel),%val(func),iosb,,,,,,,,) if (.not. status) call lib$signal(%val(status)) return end c c*****************************************************************SKIPBLOCK c subroutine skipblock (channel,nskip) ! copyright (c) 1989 by William T. Reach, University of California, Berkeley integer*4 status,sys$assign,sys$qiow,iosb(2),channel integer func*2,nskip*2 func = 38 status = sys$qiow (,%val(channel),%val(func),iosb,,, + %val(nskip),%val(1),,,,) if (.not. status) call lib$signal(%val(status)) return end c c*****************************************************************SKIPFILE c subroutine skipfile (channel,nskip) ! copyright (c) 1989 by William T. Reach, University of California, Berkeley integer*4 status,sys$assign,sys$qiow,iosb(2),channel,nskip integer func*2 func = 37 status = sys$qiow (,%val(channel),%val(func),iosb,,, + %val(nskip),%val(1),,,,) if (.not. status) call lib$signal(%val(status)) return end