PROGRAM v1access C C gfortran -o v1access v1access.f v1sub.f C C - top level program to access original URAT1 binary data zone files C - extract all items for stars in selected areas (RA, Dec box, mag range) C - formatted (ASCII) output table, various format options C - official URAT1 star designation = zone numb + run ID along zone C C use Fortran unit 13 to auxiliary data and zone files C use Fortran unit 20 for output file (or screen dump) C C 141120 NZ taken from UCAC4, adopt for URAT C 150129 NZ update for "v" and format for neg.sign item #6 IMPLICIT NONE INTEGER zmax,bmax, dimj, dima PARAMETER (zmax= 900 ! largest zone number (0.2 deg wide) . ,bmax=1440 ! numb. bins along RA for index . ,dimj= 49 ! all data items, incl. upd.pos, rec.ID . ,dima= 200000) ! max.numb. stars in a target box INTEGER alls(dima,dimj) ! all data all selected star INTEGER xi (zmax,bmax) ! index = accumul.numb.stars f(zone,RA bin) INTEGER fmto INTEGER i,j,jp,is,zn, nsr,nss,errc, uo, rah,ram,dcd,dcm . ,un,ic1, nst, ntgt, nh REAL*8 ra1,ra2,dc1,dc2,ra0,dc0,ras,dcs,mag1,mag2, newep,wra,wdc CHARACTER*40 pathz,fnxu, fnout, answer, fnlst CHARACTER*1 a1, lstfmt LOGICAL bf, was, negdec * defaults pathz = '/d01/urat/v12/' ! path to orig. UCAC4 zone files fnxu = '/d01/urat/v12/v1index.unf' fnout = 'urat1.tab' ra1 = 6.00d0 ra2 = 6.02d0 dc1 = 30.0d0 dc2 = 30.1d0 mag1 = -2.0d0 mag2 = 21.0d0 newep = 2000.0d0 wra = 0.2d0 wdc = 0.2d0 fmto = 1 was = .FALSE. ! width in arcsec, else degree bf = .FALSE. ! byte flip needed un = 13 ! Fortran unit number for zone files is = 0 ! sort by column number (0 = no sort) ic1 = 1 ! first column to sort by nst = 0 ! sum up total numb. stars output ntgt = 0 ! number of targets from list file * interactive: basic data, common for entire run of program WRITE (*,'(/a)') 'Program to access URAT1 release data' WRITE (*,'( a)') ' hit "enter" for defaults' WRITE (*,'( a)') 'option: loop through several areas possible,' WRITE (*,'( a)') ' however, output with single format to single' WRITE (*,'( a)') ' file per run of this program only' WRITE (*,'( a)') 'wrap around 24/0h RA allowed' WRITE (*,'(/a)') 'path for binary URAT1 zone files = ' WRITE (*,'(a,$)') pathz READ (*,'(a)') answer IF (answer.NE.' ') pathz = answer WRITE (*,'(/a)') 'binary data are stored in Linux-style sequence' WRITE (*,'( a)') 'some computers (PC-style) need a byteflip' WRITE (*,'(a,l1,a,$)') 'byte flip (true/false) ? ',bf,' ' READ (*,'(a)') a1 IF (a1.NE.' ') READ (a1,*) bf WRITE (*,'(/a)') 'full path and file name of unform. index =' WRITE (*,'(a,$)') fnxu READ (*,'(a)') answer IF (answer.NE.' ') fnxu = answer WRITE (*,'(/a)') ' fmt = 1 = all original integers' WRITE (*,'( a)') ' fmt = 2 = updated RA, Dec (degree) short' WRITE (*,'( a)') ' fmt = 3 = updated RA, Dec (hms) photom..' WRITE (*,' (a)') ' fmt = ... user defined in subr. put_stars' WRITE (*,'(a,i2,a,$)') 'format for output = ',fmto,' ' READ (*,'(a)') answer IF (answer.NE.' ') READ (answer,*) fmto WRITE (*,'(/a)') 'file name output table or "s" = screen dump ' WRITE (*,'(a,$)') fnout READ (*,'(a)') answer IF (answer.NE.' ') fnout = answer WRITE (*,'(/a)') . 'select magnitude range (URAT mags, -2..21 = all):' WRITE (*,'(2f7.2,a,$)') mag1, mag2,' ' READ (*,'(a)') answer IF (answer.NE.' ') READ (answer,*) mag1,mag2 IF (fmto.GT.1) THEN WRITE (*,'(/a)') . 'update positions for proper motions to new epoch' WRITE (*,'( a)') '2013.5 is default = mean URAT1 obs.epoch' WRITE (*,'(f7.2,a,$)') newep,' ' READ (*,'(a)') answer IF (answer.NE.' ') READ (answer,*) newep ENDIF WRITE (*,'(/a)') 'see readme.txt for all items' WRITE (*,'( a)') ' 0 = no sort = default' WRITE (*,'( a)') ' 1 = sort by col. 1 = RA ' WRITE (*,'( a)') ' 2 = sort by col. 2 = Dec' WRITE (*,'( a)') ' 8 = sort by URAT mag' WRITE (*,'( a)') ' ...' WRITE (*,'( a)') ' 49 = sort by URAT ID number' WRITE (*,'(a,i2,a,$)') 'sort by column number: ',is,' ' READ (*,'(a)') answer IF (answer.NE.' ') READ (answer,*) is * read index file (900 x 1440 x 4 byte) WRITE (*,'(/a,a)') 'open index file = ',fnxu OPEN (un,FILE=fnxu,ACCESS='direct',RECL=5184000) READ (un,REC=1) xi CLOSE(un) WRITE (*,'(a)') 'index read successfully, sample data:' DO zn=600,604 WRITE (*,'(a,6i6)') 'zn, xi...=',zn,(xi(zn,j),j=1,5) ENDDO IF (bf) THEN CALL xi_byte_flip (xi,zmax,bmax) ! apply byte flip if needed WRITE (*,'(a)') 'sample of index data after byte flip:' DO zn=600,604 WRITE (*,'(a,6i6)') 'zn, xi...=',zn,(xi(zn,j),j=1,5) ENDDO ENDIF * prepare table output IF (fnout.EQ.'s') THEN uo = 6 ELSE uo = 20 OPEN (uo,FILE=fnout) ENDIF * input list or interactive WRITE (*,'(/a,$)') . 'name of list with centers or blank for interactive: ' READ (*,'(a)') fnlst WRITE (*,'(a,$)') '(h) = hms format, (d) = decimal degree ' READ (*,'(a)') lstfmt IF (fnlst.NE.' ') THEN OPEN (12,FILE=fnlst) WRITE (*,'(/a,2f7.3,a,$)') . 'width of box RA*cosD, Dec [degree] =',wra,wdc,' ' READ (*,'(a)') answer IF (answer.NE.' ') READ (answer,*) wra,wdc ENDIF * loop boxes 101 CONTINUE IF (fnlst.EQ.' ') THEN WRITE (*,'(/a,$)') '== new box: r=range, c=center, q=quit ' READ (*,'(a)') a1 IF (a1.EQ.'r'.OR.a1.EQ.'R') THEN WRITE (*,'(a,2f8.4,a,$)') 'RA (hour) = ',ra1,ra2,' ' READ (*,'(a)') answer IF (answer.NE.' ') READ (answer,*) ra1,ra2 WRITE (*,'(a,2f8.3,a,$)') 'Dec (deg) = ',dc1,dc2,' ' READ (*,'(a)') answer IF (answer.NE.' ') READ (answer,*) dc1,dc2 ELSEIF (a1.EQ.'c'.OR.a1.EQ.'C') THEN WRITE (*,'(2a)') 'width in degree or enter "a" for arcsec' . ,' after values' WRITE (*,'(a,2f7.3,a,$)') . 'width of box RA*cosD, Dec =',wra,wdc,' ' READ (*,'(a)') answer j = INDEX(answer,'a') IF (j.GT.0) THEN was = .TRUE. answer (j:j) = ' ' ELSE was = .FALSE. ENDIF IF (answer.NE.' ') THEN READ (answer,*) wra,wdc IF (was) THEN wra = wra / 3.6d3 ! arcsec to degree wdc = wdc / 3.6d3 ENDIF ENDIF IF (wra.LE.0.0d0.OR.wdc.LE.0.0d0) GOTO 101 IF (lstfmt.EQ.'h') THEN WRITE (*,'(a,$)') 'center RA (hms), Dec (dms) = ' READ (*,'(a)') answer CALL str2deg (answer,ra0,dc0) CALL calc_box (ra0,dc0,wra,wdc,ra1,ra2,dc1,dc2) WRITE (*,'(2f10.6,2f10.5)') ra1,ra2,dc1,dc2 ELSE WRITE (*,'(a,$)') 'center RA (deg), Dec (deg) = ' READ (*,*) ra0,dc0 ! degree WRITE (*,'(a,4f12.6)') 'ra0,dc0 = ',ra0,dc0, wra,wdc CALL calc_box (ra0,dc0,wra,wdc,ra1,ra2,dc1,dc2) ENDIF ELSE GOTO 99 ! exit interactive loop ENDIF ELSE ! read centers from list file IF (lstfmt.EQ.'h') THEN READ (12,'(a)',END=99) answer ! hms format CALL str2deg (answer,ra0,dc0) CALL calc_box (ra0,dc0,wra,wdc,ra1,ra2,dc1,dc2) ELSE READ (12,*,END=99) ra0,dc0 ! degree WRITE (*,'(a,4f12.6)') 'ra0,dc0 = ',ra0,dc0, wra,wdc CALL calc_box (ra0,dc0,wra,wdc,ra1,ra2,dc1,dc2) ENDIF ntgt = ntgt + 1 WRITE (*,'(/a,2f10.6,2f10.5)') 'new field:', ra1,ra2,dc1,dc2 ENDIF * output info IF (fnlst.EQ.' ') THEN WRITE (uo,'(/a, f10.3)') 'epoch = ',newep WRITE (uo,'( a, i10 )') 'sort item = ',is WRITE (uo,'( a, i10 )') 'output fmt = ',fmto WRITE (uo,'( a,2f10.2)') 'mag range = ',mag1,mag2 WRITE (uo,'( a,2f10.2)') 'box width = ',wra, wdc WRITE (uo,'( a,f10.6,f10.5)') 'min RA, Dec= ',ra1,dc1 WRITE (uo,'( a,f10.6,f10.5)') 'max RA, Dec= ',ra2,dc2 ENDIF * get stars CALL get_stars (bf,ra1,ra2,dc1,dc2,mag1,mag2 . ,newep,pathz,xi,zmax,bmax,dima,dimj . ,alls,nss,nsr,errc) WRITE (*,'(a,3i7)') . 'stars read, selected, eofile = ',nsr,nss,errc * sort IF (is.GE.1.AND.is.LE.dimj) THEN WRITE (*,'(a,i6)') 'start sort by column = ',is CALL sorti (alls,is,dimj,ic1,nss,dima,dimj) ENDIF * output CALL put_stars (ntgt,newep,alls,dima,dimj,nss,fmto,uo) WRITE (*,'(a,i6,a)') 'output complete for ',nss,' stars' nst = nst + nss GOTO 101 ! more boxes * the end 901 WRITE (*,'(a)') 'error in reading index file' STOP 99 CONTINUE WRITE (*,'(/a,i8)') 'total numb. stars output = ',nst WRITE (*,'(a/)') 'thanks for using , have a nice day!' END ! main ************************************************************************ SUBROUTINE str2deg (answer,ra0,dc0) C C input : answer = string with RA, Dec in hms format C output: ra0,dc0= RA, Dec [degree] IMPLICIT NONE CHARACTER*(*) answer REAL*8 ra0,dc0, ras,dcs INTEGER rah,ram,dcd,dcm, j LOGICAL negdec j = INDEX (answer,'-') IF (j.GT.0) THEN negdec = .TRUE. answer (j:j) = ' ' ELSE negdec = .FALSE. ENDIF READ (answer,*) rah,ram,ras, dcd,dcm,dcs ra0 = rah + ram/60.0d0 + ras/3.6d3 dc0 = dcd + dcm/60.0d0 + dcs/3.6d3 IF (negdec) dc0 = -dc0 IF (ra0.LT. 0.0d0.OR.ra0.GT.360.0d0.OR. . dc0.LT.-90.0d0.OR.dc0.GT. 90.0d0) THEN WRITE (*,'(/a)') 'error ' WRITE (*,'(a)') answer WRITE (*,'(2f12.6)') ra0,dc0 STOP ENDIF END ! subr. ************************************************************************ SUBROUTINE calc_box (rad,dcd,wra,wdc,ra1,ra2,dc1,dc2) C C input : rad,dcd= RA, Dec [degree] C wra = width of box [degree] along RA (*cosDec) C wdc = width of box [degree] along Dec C output: ra1,ra2= range in RA [hour] C dc1,dc2= range in Declination [degree] IMPLICIT NONE REAL*8 rad,dcd, wra,wdc, ra1,ra2,dc1,dc2, rah REAL*8 degrad, cosdec, hwr,hwd, ras,dcs degrad = DATAN (1.0d0) / (4.5d1) ! degree to radian cosdec = DCOS (dcd * degrad) rah = rad / 15.0d0 ! degree -> hour hwr = wra / (cosdec * 30.0d0) ! half width in hour for RA hwd = wdc / 2.0d0 ! half width in Dec ra1 = rah - hwr ra2 = rah + hwr IF (ra1.LT. 0.0d0) ra1 = ra1 + 24.0d0 IF (ra2.GT.24.0d0) ra2 = ra2 - 24.0d0 dc1 = dcd - hwd dc2 = dcd + hwd END ! subr. ************************************************************************ SUBROUTINE get_stars (bf,ra1,ra2,dc1,dc2,mag1,mag2 . ,newep,pathz,xi,zmax,bmax,dima,dimj . ,alls,nss,nsr,errc) C C get all data items for all stars within specified RA,Dec box and mag range C C input : bf = .TRUE. if flip of bytes is required C ra1,ra2 = RA range (hour, decimal) C dc1,dc2 = declination range (degree, decimal) C mag1,mag2= range URAT magnitudes to pick stars C newep = requested epoch for positions (year) C pathz = path name for URAT1 zone files C xi (zmax,bmax) = index array = f(zone, RA bin) C index runs along indiv. zone only = last star ID of bin C zmax = largest zone number (0.2 deg steps from South Pole) C bmax = largest bin along RA (0.25 deg steps = 1 arcmin) C dima = max. number of stars to retrieve C dimj = number of items per star incl. star ID number C C output: alls(dima,dimj) = all items for all retrieved stars C nss = number of stars selected = within specified box C nsr = number of star records read from zone files C errc = number of errors while reading zone files IMPLICIT NONE INTEGER dimc PARAMETER (dimc = 45) ! numb. columns (data items from file) LOGICAL bf REAL*8 ra1,ra2,dc1,dc2, mag1,mag2, newep CHARACTER*(*) pathz INTEGER zmax,bmax,dima,dimj, nb, nsr,nss,errc INTEGER xi (zmax,bmax) ! accumul.numb.stars = f(zone,RA bin) INTEGER alls(dima,dimj) ! all integer data each star INTEGER idat(dimc) ! data items individual star INTEGER i,j,zn,rr,uz, d1m,d2m,z1,z2,nz, mi1,mi2 . ,r1m(2),r2m(2),i1(2),i2(2),nr, i0,j1,j2, recn . ,ran,spdn, spn, sd1,sd2 REAL*8 masrad CHARACTER*40 fnzone LOGICAL errflg * prepare IF (dimj.LT.49) THEN WRITE (*,'(a)') 'error need dimj >= 49' STOP ENDIF mi1 = IDNINT (mag1 * 1.0d3) ! from R*8 magnitudes mi2 = IDNINT (mag2 * 1.0d3) ! to I*4 [milli-mag] masrad = DATAN (1.0d0) / (4.5d1*3.6d6) ! mas to radian * calculate range CALL get_zone_range (dc1,dc2,zmax,d1m,d2m,z1,z2,nz) WRITE (*,'(a,2i5,2i11)') '=== z1,2,d1,2 = ',z1,z2,d1m,d2m sd1 = d1m + 324000000 ! Dec -> SPD [mas] sd2 = d2m + 324000000 IF (nz.LT.1) THEN nss = 0 RETURN ! exit, declination invalid ENDIF CALL get_ra_range (ra1,ra2,r1m,r2m,i1,i2,nr) WRITE (*,'(a,4i11)') '=== r1m,2 = ',r1m,r2m WRITE (*,'(a,4i11)') '=== i1,i2 = ',i1,i2 * initialize uz = 13 ! unit number for zone files nsr = 0 nss = 0 errc= 0 * loop all zones DO zn = z1,z2 WRITE (*,'(a,i3)') 'open file for zone = ',zn CALL open_zfile (pathz,uz,zn,fnzone) WRITE (*,'(a)') fnzone * rr = 1 or 2 = numb. of RA ranges, normal=1, else crossover 24/0 hour * i1,i2 = 1,...,1440 = bin number (0.25 deg) along RA * j1,j2 = 1,... many = record number on individ. zone file (1 record = 1 star) * xi(..)= largest record number of this RA bin or previous non-zero bin DO rr = 1,nr ! 1 or 2 RA ranges possible errflg = .FALSE. ! reset, no end of file encountered yet i0 = i1(rr) - 1 ! previous RA bin IF (i0.LT.1) THEN ! data request begins at begin 1st RA bin j1 = 1 ! = very first possible star on zone file ELSE ! normal case j1 = xi(zn,i0) + 1 ! 1st star to use = prev.bin last + 1 ENDIF ! (xi = contains last rec.# in bin) j2 = xi(zn,i2(rr)) ! last rec.# = end of last bin to use WRITE (*,'(a,i3,4i5,2i7)') . 'zn, rr,i0,i1,i2, j1,j2 = ',zn,rr,i0,i1(rr),i2(rr),j1,j2 DO recn= j1,j2 CALL getistar (uz,recn,bf,errflg,idat,dimc) IF (errflg) THEN errc = errc + 1 ! count read errors = end of files encountered GOTO 91 ELSE nsr = nsr + 1 ! count number of stars read IF (idat(8).GE.mi1.AND.idat(8).LE.mi2) THEN ! URAT mag range IF (idat(1).GE.r1m(rr).AND.idat(1).LE.r2m(rr).AND. . idat(2).GE.sd1.AND.idat(2).LE.sd2) THEN nss = nss + 1 IF (nss.GT.dima) THEN WRITE (*,'(/a,i8)') . 'WARNING: too many stars, take first dima =',dima RETURN ENDIF DO j=1,dimc ! all data columns alls (nss,j) = idat(j) ENDDO CALL apply_pm (idat,dimc,newep,masrad, ran,spdn,spn) alls (nss,46) = ran alls (nss,47) = spdn alls (nss,48) = spn ! = same error for both coord. alls (nss,49) = zn * 1000000 + recn ENDIF ! within box ENDIF ! within mag range ENDIF ! case read o.k. or error ENDDO ! all stars within a zone 91 CONTINUE ! end of zone file WRITE (*,'(a,i2,2i6)') '=== rr,nsr,nss = ',rr,nsr,nss ENDDO ! 1 or 2 RA ranges ENDDO ! all zones END ! subr. ************************************************************************ SUBROUTINE put_stars (ntgt,newep,alls,dima,dimj,nss,fmto,uo) C C input : ntgt = target number (from list) or 0 C newep = target epoch of position data C alls = array with data of all selected stars, all columns C dima = max.numb. of stars allowed in target area C dimj = number of items per star C nss = number of selected stars C fmto = format modus for output (= 1 to 3 here) C uo = Fortran unit number for output file C output : items to file with unit number uo (can be screen = std.out) C header for ntgt <= 1 IMPLICIT NONE INTEGER dima,dimj,nss, fmto,uo INTEGER alls(dima,dimj) INTEGER ntgt,i,j, nrj,spd,ep97, nmag,area,wdsf, dnn,hstf,badf REAL*8 newep,magm, epr,epd, pmr,pmrc,pmd . ,ran,dcn, magj, magh, magk, epmr,epmd . ,magb,magv,magg,magr,magi CHARACTER*13 cra,cdc * check IF (fmto.LT.1.OR.fmto.GT.3) THEN WRITE (*,'(a,i3)') '** invalid format fmto = ',fmto RETURN ENDIF * header line IF (ntgt.LE.0) THEN ! header for no-list input IF (fmto.EQ.1) THEN WRITE (uo,'(/2i10,2i4,i3,i4,i6,i6,i4,i3,i2,4i4,2i6,i4 . ,2i3,i11,3i6,3i5,6i2,5i6,5i5,2i4,i11,i10,i5,i10)') . (j,j=1,49) WRITE (uo,'(6a)') ' RA SPD' . ,' ss sm ns nu epoc mag err nm r not nou got gou' . ,' pmr pmd pme m2 ma 2MASS ID J H K' . ,' eJ eH eK c c c q q q B V g r' . ,' i eB eV eg er ei nn no' . ,' RAn SPDn sx zn rec' ELSEIF (fmto.EQ.2) THEN WRITE (uo,'(/2a)') ' RA (deg) DE (deg)' . ,' err mag o.epoch PMrac PMdec ePM zn rec' ELSEIF (fmto.EQ.3) THEN WRITE (uo,'(/3a)') ' RA DE' . ,' err nu Umag Jmag Hmag Kmag' . , ' Bmag Vmag gmag rmag imag zn rec' ENDIF ELSEIF (ntgt.EQ.1) THEN ! header for list input IF (fmto.EQ.1) THEN WRITE (uo,'(/i4,i11,i10,2i4,i3,i4,i6,i6,i4,i3,i2,4i4,2i6,i4 . ,2i3,i11,3i6,3i5,6i2,5i6,5i5,2i4,i11,i10,i5,i10)') . (j,j=0,49) WRITE (uo,'(6a)') ' nt RA SPD' . ,' ss sm ns nu epoc mag err nm r not nou got gou' . ,' pmr pmd pme m2 ma 2MASS ID J H K' . ,' eJ eH eK c c c q q q B V g r' . ,' i eB eV eg er ei nn no' . ,' RAn SPDn sx zn rec' ELSEIF (fmto.EQ.2) THEN WRITE (uo,'(/2a)') ' nr RA (deg) DE (deg)' . ,' err mag o.epoch PMrac PMdec ePM zn rec' ELSEIF (fmto.EQ.3) THEN WRITE (uo,'(/3a)') ' nt RA DE' . ,' err nu Umag Jmag Hmag Kmag' . , ' Bmag Vmag gmag rmag imag zn rec' ENDIF ENDIF ! header only for 1st target from list or interactive cases * table with data DO i=1,nss IF (fmto.EQ.1) THEN ! all original integer IF (ntgt.LE.0) THEN WRITE (uo,'(2i10,2i4,i3,i4,i6,i6,i4,i3,i2,4i4,2i6,i4 . ,2i3,i11,3i6,3i5,6i2,5i6,5i5,2i4,i11,i10,i5,i10)') . (alls(i,j),j=1,49) ELSE WRITE (uo,'(i4,i11,i10,2i4,i3,i4,i6,i6,i4,i3,i2,4i4,2i6,i4 . ,2i3,i11,3i6,3i5,6i2,5i6,5i5,2i4,i11,i10,i5,i10)') . ntgt,(alls(i,j),j=1,49) ENDIF ELSEIF (fmto.EQ.2) THEN ran = alls(i,46) / 3.6d6 ! updated mas -> degree dcn = alls(i,47) / 3.6d6 - 90.0d0 ! updated mas -> degree magm= alls(i, 8) * 1.0d-3 epr = alls(i, 7) * 1.0d-3 + 2000.0d0 pmrc= alls(i,16) * 0.1d0 pmd = alls(i,17) * 0.1d0 epmr= alls(i,18) * 0.1d0 IF (ntgt.LE.0) THEN WRITE (uo,'(f11.7,f12.7,i4, f7.3,f9.3,2f7.1,f5.1,i10)') . ,ran,dcn,alls(i,48), magm,epr,pmrc,pmd,epmr, alls(i,49) ELSE WRITE (uo,'(i4,2f12.7,i4, f7.3,f9.3,2f7.1,f5.1,i10)') . ntgt,ran,dcn,alls(i,48), magm,epr,pmrc,pmd,epmr, alls(i,49) ENDIF ELSEIF (fmto.EQ.3) THEN ran = alls(i,46) * 1.0d-3 ! updated pos. mas -> arcsec dcn = alls(i,47) * 1.0d-3 - 324.0d3 ! updated pos. mas -> arcsec magm= alls(i, 8) * 1.0d-3 magj= alls(i,22) * 1.0d-3 ! from 2MASS magh= alls(i,23) * 1.0d-3 magk= alls(i,24) * 1.0d-3 magb= alls(i,34) * 1.0d-3 ! APASS mags magv= alls(i,35) * 1.0d-3 magg= alls(i,36) * 1.0d-3 magr= alls(i,37) * 1.0d-3 magi= alls(i,38) * 1.0d-3 CALL as2hms (ran,dcn,cra,cdc) IF (ntgt.LE.0) THEN WRITE (uo,'(a,1x,a,2i4,9f7.3,i10)') . cra,cdc,alls(i,48),alls(i,13),magm . ,magj,magh,magk,magb,magv,magg,magr,magi, alls(i,49) ELSE WRITE (uo,'(i4,1x,a,1x,a,2i4,9f7.3,i10)') . ntgt,cra,cdc,alls(i,48),alls(i,13),magm . ,magj,magh,magk,magb,magv,magg,magr,magi, alls(i,49) ENDIF ENDIF ! case fmto ENDDO ! all stars END ! subr. ************************************************************************ SUBROUTINE pos_up (ra,spd,pmra,pmdc,oldep,newep, ran,spdn) C C update position for proper motion with epoch difference C C input : ra, spd = RA, SPD (mas) at input epoch C pmra,pmdc = proper motion RA, Dec (0.1 mas/yr) C no cos(Dec) for pmra, i.e. large values near pole C oldep = old epoch (year) of input data C newep = new epoch (year) C output: ran,spdn = RA, SPD (mas) at new epoch IMPLICIT NONE INTEGER ra,spd,pmra,pmdc,ran,spdn REAL*8 oldep, newep, dt, dra, ddc dt = newep - oldep dra = DFLOAT(pmra) * 0.1d0 * dt ! PM RA (without cos Dec) ddc = DFLOAT(pmdc) * 0.1d0 * dt C WRITE (33,'(/a,3f11.3)') 'dt,dra,ddc = ',dt,dra,ddc ran = ra + IDNINT (dra) IF (ran.GT.1296000000) ran = ran - 1296000000 ! 24 hour in mas IF (ran.LT. 0) ran = ran + 1296000000 ! restrict to 0..24 range spdn= spd + IDNINT (ddc) ! assume no jump over pole END ! subr. ************************************************************************ SUBROUTINE pos_error (sigc,spm, cep,newep, spn) C C calculate position error at new epoch per coordinate C assume same error for RA, Dec coordinates (URAT1 data) C if no PM given (PM error = 90 mas/yr) position error is set to 900 mas C C input : sigc = position error (mas) at central epoch C spm = error of proper motion component (0.1 mas/yr) C cep = central epoch (1/100 yr after 1900) C newep= requested epoch for position error (year) C output: spn = position error at new epoch (mas) IMPLICIT NONE INTEGER sigc, spm, spn REAL*8 cep, newep, dt dt = newep - cep ! year IF (spm.GE.900) THEN ! no PM known spn = 900 ! default position error at new epoch ELSE ! normal case, use valid PM error spn = IDNINT (DSQRT (sigc**2 + (spm * 0.1d0 * dt)**2)) ENDIF END ! subr. ************************************************************************ SUBROUTINE apply_pm (idat,dimc,newep,masrad, ran,spdn,spn) C C input : idat = integer vector with original URAT1 data for 1 star C dimc = dimension of idata vector C newep= requested new epoch (year) C masrad = convert mas to radian C output: ran,spdn = new position, updated for proper motion (mas) C spn = position error at new epoch, per coord. (mas) IMPLICIT NONE INTEGER dimc, idat(dimc), ran,spdn,spn, pmra REAL*8 oldep, newep, masrad, dec,cosd oldep= idat(7)*1.0d-3 + 2000.0d0 ! orig. URAT obs. epoch dec = DFLOAT(idat(2)-324000000) ! SPD [mas] -> Dec [mas] cosd = DCOS (dec*masrad) pmra = IDNINT (idat(16)/cosd) ! PM RA (without cosDec factor) CALL pos_up (idat(1),idat(2),pmra,idat(17),oldep,newep,ran,spdn) CALL pos_error (idat(3),idat(18), oldep,newep, spn) C WRITE (33,'(a,2i11,2i6)') 'ra,spd,err,mag = ' C . ,idat(1),idat(2),idat(3),idat(8) C WRITE (33,'(a,2i11,2i6)') 'pmr,pmd,pme,epo= ' C . ,idat(16),idat(17),idat(18),idat(7) C WRITE (33,'(a,2f11.3,f11.0,f9.4,i8)') 'oldep,newep,dec= ' C . ,oldep,newep,dec,cosd,pmra C WRITE (33,'(a,2i11,2i6)') 'new ra,spd,err = ' C . ,ran,spdn,spn END ! subr. ************************************************************************