!----------------------------------------------------------------------------------------------------------------------------------! ! PROGRAM PPMXL ! ! ! ! Creates a direct-access version of the PPMXL catalog, splitting the catalog into one-degree-wide declination zones sorted by ! ! right ascension and with bad magnitudes suppressed ! ! ! ! Copyright (C) 2011 by David J. Tholen ! ! Last update: 2011 Dec 17 ! !----------------------------------------------------------------------------------------------------------------------------------! PROGRAM PPMXL IMPLICIT NONE REAL (KIND=4) :: B1Mag ! first epoch B magnitude REAL (KIND=4) :: B2Mag ! second epoch B magnitude REAL (KIND=4) :: BDiff ! difference in B magnitudes REAL (KIND=4) :: BMag ! average B magnitude INTEGER :: BMagArray (-110:670) ! blue magnitude array for histogram REAL (KIND=8) :: Dec ! declination in degrees INTEGER, PARAMETER :: DegToMAS = 3600000 ! degrees to milliarcsec conversion factor REAL (KIND=4) :: EpochDec ! declination mean epoch in year REAL (KIND=4) :: EpochRA ! right ascension mean epoch in year INTEGER :: ErrorCode ! error code REAL (KIND=4) :: EWErr ! east-west error in degrees (E.-W. = R.A. times cosine DEC.) REAL (KIND=4) :: EWPM ! east-west proper motion in degrees per year REAL (KIND=4) :: EWPMErr ! east-west proper motion error in degrees per year CHARACTER (LEN= 12) :: Field ! twelve character subset of Line CHARACTER (LEN= 7) :: File ! filename INTEGER :: Flags ! flags INTEGER :: I ! character substring start index or DO loop index INTEGER (KIND=2) :: IBMag ! integer version of chosen B magnitude in hundredths of a magnitude INTEGER :: IDec ! integer version of declination in milliarcsec INTEGER (KIND=2) :: IEpochDec ! integer declination epoch in hundredths of a year from 1900 INTEGER (KIND=2) :: IEpochRA ! integer right ascension epoch in hundredths of a year from 1900 INTEGER (KIND=1) :: IEWErr ! east-west error in milliarcsec minus 128 INTEGER :: IEWPM ! east-west proper motion in 0.1 milliarcsec per year INTEGER (KIND=2) :: IEWPMErr ! east-west proper motion error in 0.1 milliarcsec per year INTEGER (KIND=1) :: INSErr ! north-south error in milliarcsec minus 128 INTEGER :: INSPM ! north-south proper motion in 0.1 milliarcsec per year INTEGER (KIND=2) :: INSPMErr ! north-south proper motion error in 0.1 milliarcsec per year INTEGER :: IRA ! integer version of right ascension in milliarcsec INTEGER (KIND=2) :: IRMag ! integer version of chosen R magnitude in hundredths of a magnitude INTEGER :: J ! character substring end index INTEGER :: K ! array index INTEGER :: Limit ! threshold for generating another entry in the accelerator file CHARACTER (LEN=200) :: Line ! one PPMXL catalog record REAL (KIND=4) :: MaxEpoch = 0.0 ! maximum mean epoch REAL (KIND=4) :: MaxErr = 0.0 ! maximum positional error REAL (KIND=4) :: MaxMag = 0.0 ! faintest magnitude REAL (KIND=4) :: MaxPM = 0.0 ! maximum proper motion REAL (KIND=4) :: MaxPMErr = 0.0 ! maximum proper motion error REAL (KIND=4) :: MinEpoch = 2011.0 ! minimum mean epoch REAL (KIND=4) :: MinErr = 9.9 ! minimum positional error REAL (KIND=4) :: MinMag = 9.9 ! brightest magnitude REAL (KIND=4) :: MinPM = 9.9 ! minimum proper motion REAL (KIND=4) :: MinPMErr = 9.9 ! minimum proper motion error CHARACTER (LEN= 4), PARAMETER :: None = 'None' ! the string "None" REAL (KIND=4) :: NSErr ! north-south error in degrees REAL (KIND=4) :: NSPM ! north-south proper motion in degrees per year REAL (KIND=4) :: NSPMErr ! north-south proper motion error in degrees per year INTEGER :: NumEpochD ! number of entries with epoch differences INTEGER :: NumHighPM ! number of high proper motion entries INTEGER (KIND=8) :: NumMags ! total number of magnitude entries (non-None, both B and R) INTEGER :: NumNoMags ! number of sources without any magnitude information INTEGER :: NumVar ! number of variable sources CHARACTER (LEN= 1), PARAMETER :: Pipe = '|' ! the pipe character INTEGER :: Progress = 0 ! progress counter INTEGER, ALLOCATABLE :: Ptr (:) ! pointer array REAL (KIND=4) :: R1Mag ! first epoch R magnitude REAL (KIND=4) :: R2Mag ! second epoch R magnitude REAL (KIND=8) :: RA ! right ascension in degrees INTEGER, ALLOCATABLE :: RAArray (:) ! array of right ascensions in a zone REAL (KIND=4) :: RDiff ! difference in R magnitudes CHARACTER (LEN=63), ALLOCATABLE :: Rest (:) ! string to hold the rest of a record (everything but R.A.) REAL (KIND=4) :: RMag ! average R magnitude INTEGER :: RMagArray (-110:670) ! red magnitude array for histogram INTEGER :: Size ! number of entries in a zone INTEGER :: Sum ! sum of zone counts (should match total number of records) LOGICAL :: Variable ! variability flag INTEGER :: ZoneCount (0:179) ! number of entries in each declination zone INTEGER :: ZoneNum ! zone number !----------------------------------------------------------------------------------------------------------------------------------! ! Open input monolithic PPMXL catalog file and ancillary output files ! ! ! ! Do not use logical unit numbers 0 to 179, which are used by the 180 declination zones ! ! The anomalies file is used to save unusual catalog entries ! ! The stats file is used to report various catalog statistics ! ! The histogram.dat file is used for the magnitude distributions ! !----------------------------------------------------------------------------------------------------------------------------------! OPEN (199,FILE='ppmxl') OPEN (198,FILE='anomalies') OPEN (197,FILE='stats') OPEN (196,FILE='histogram.dat') !----------------------------------------------------------------------------------------------------------------------------------! ! Initialize the magnitude arrays and accumulators ! !----------------------------------------------------------------------------------------------------------------------------------! BMagArray = 0 RMagArray = 0 NumEpochD = 0 NumHighPM = 0 NumMags = 0 NumNoMags = 0 NumVar = 0 !----------------------------------------------------------------------------------------------------------------------------------! ! Open raw output file for all 180 declination zones and initialize number of entries in each zone ! !----------------------------------------------------------------------------------------------------------------------------------! DO I=0,179,1 WRITE (File,'(I3.3,A4)') I,'.raw' OPEN (I,FILE=File) ZoneCount(I) = 0 END DO !----------------------------------------------------------------------------------------------------------------------------------! ! Begin main loop to read all records from the monolithic PPMLX file ! !----------------------------------------------------------------------------------------------------------------------------------! DO READ (199,'(A)',IOSTAT=ErrorCode) Line IF (ErrorCode < 0) EXIT I = INDEX(Line(1:),Pipe) + 1 ! find lower limit of right ascension field J = INDEX(Line(I:),Pipe) + I - 2 ! find upper limit of right ascension field Field = Line(I:J) READ (Field,'(F12.0)') RA ! extract right ascenion I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 Field = Line(I:J) READ (Field,'(F12.0)') Dec ! extract declination I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 Field = Line(I:J) READ (Field,'(F12.0)') EWErr ! extract east-west error (R.A. error times cosine DEC.) I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 Field = Line(I:J) READ (Field,'(F12.0)') NSErr ! extract north-south error (DEC. error) I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 Field = Line(I:J) READ (Field,'(F12.0)') EWPM ! extract east-west proper motion I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 Field = Line(I:J) READ (Field,'(F12.0)') NSPM ! extract north-south proper motion I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 Field = Line(I:J) READ (Field,'(F12.0)') EWPMErr ! extract east-west proper motion error I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 Field = Line(I:J) READ (Field,'(F12.0)') NSPMErr ! extract north-south proper motion error I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 ! skip number of observations I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 ! extract RA epoch Field = Line(I:J) READ (Field,'(F7.0)') EpochRA I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 ! extract Dec epoch Field = Line(I:J) READ (Field,'(F7.0)') EpochDec I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 ! skip J mag I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 ! skip J mag err I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 ! skip H mag I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 ! skip H mag err I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 ! skip K mag I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 ! skip K mag err I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 Field = Line(I:J) IF (Field .NE. None) THEN ! if first epoch B magnitude is present READ (Field,'(F12.0)') B1Mag ! extract first epoch B magnitude IF (B1Mag > MaxMag) MaxMag = B1Mag ! does it hold the faint record for B? IF (B1Mag < MinMag) MinMag = B1Mag ! does it hold the bright record for B? ! Save anomalously bright or faint magnitudes to the anomalies file IF (B1Mag > 24.5) WRITE (198,'(2A,F6.2)') Line(1:35),' B1 ',B1Mag IF (B1Mag < -2.0) WRITE (198,'(2A,F6.2)') Line(1:35),' B1 ',B1Mag K = INT(10*B1Mag) BMagArray(K) = BMagArray(K) + 1 ! add to blue magnitude histogram array NumMags = NumMags + 1 ! increment number of magnitudes ELSE B1Mag = 99.99 ! if no first epoch B magnitude is present, assign a large default END IF I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 Field = Line(I:J) IF (Field .NE. None) THEN ! if second epoch B magnitude is present READ (Field,'(F12.0)') B2Mag ! extract second epoch B magnitude IF (B2Mag > MaxMag) MaxMag = B2Mag ! does it hold the faint record for B? IF (B2Mag < MinMag) MinMag = B2Mag ! does it hold the bright record for B? ! Save anomalously bright or faint magnitudes to the anomalies file IF (B2Mag > 24.5) WRITE (198,'(2A,F6.2)') Line(1:35),' B2 ',B2Mag IF (B2Mag < -2.0) WRITE (198,'(2A,F6.2)') Line(1:35),' B2 ',B2Mag K = INT(10*B2Mag) BMagArray(K) = BMagArray(K) + 1 ! add to blue magnitude histogram array NumMags = NumMags + 1 ! increment number of magnitudes ELSE B2Mag = 99.99 ! if not second epoch B magnitude is present, assign a large default END IF I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 Field = Line(I:J) IF (Field .NE. None) THEN ! if first epoch R magnitude is present READ (Field,'(F12.0)') R1Mag ! extract first epoch R magnitude IF (R1Mag > MaxMag) MaxMag = R1Mag ! does it hold the faint record for R? IF (R1Mag < MinMag) MinMag = R1Mag ! does it hold the bright record for R? ! Save anomalously bright or faint magnitudes to the anomalies file IF (R1Mag > 23.5) WRITE (198,'(2A,F6.2)') Line(1:35),' R1 ',R1Mag IF (R1Mag < -2.0) WRITE (198,'(2A,F6.2)') Line(1:35),' R1 ',R1Mag K = INT(10*R1Mag) RMagArray(K) = RMagArray(K) + 1 ! add to red magnitude histogram array NumMags = NumMags + 1 ! increment number of magnitudes ELSE R1Mag = 99.99 ! if no first epoch R magnitude is present, assign a large default END IF I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 Field = Line(I:J) IF (Field .NE. None) THEN ! if second epoch R magnitude is present READ (Field,'(F12.0)') R2Mag ! extract second epoch R magnitude IF (R2Mag > MaxMag) MaxMag = R2Mag ! does it hold the faint record for R? IF (R2Mag < MinMag) MinMag = R2Mag ! does it hold the bright record for R? ! Save anomalously bright or faint magnitudes to the anomalies file IF (R2Mag > 23.5) WRITE (198,'(2A,F6.2)') Line(1:35),' R2 ',R2Mag IF (R2Mag < -2.0) WRITE (198,'(2A,F6.2)') Line(1:35),' R2 ',R2Mag K = INT(10*R2Mag) RMagArray(K) = RMagArray(K) + 1 ! add to red magnitude histogram array NumMags = NumMags + 1 ! increment number of magnitudes ELSE R2Mag = 99.99 ! if no second epoch R magnitude is present, assign a large default END IF I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 ! skip field I = J + 2 J = INDEX(Line(I:),Pipe) + I - 2 ! skip field I = J + 2 Field = Line(I:) READ (Field,'(I10)') Flags ! extract flags !----------------------------------------------------------------------------------------------------------------------------------! ! Finished reading record and extracting useful quantities, now convert values ! !----------------------------------------------------------------------------------------------------------------------------------! ZoneNum = Dec + 90.0_8 ! determine declination zone (truncates real to integer; use one degree wide zones) ! Find largest and smallest values of position error, proper motion, and proper motion error ! so that an encoding algorithm can be chosen IF (EWErr > MaxErr) MaxErr = EWErr IF (NSErr > MaxErr) MaxErr = NSErr IF (EWErr < MinErr) MinErr = EWErr IF (NSErr < MinErr) MinErr = NSErr IF (EWPM > MaxPM) MaxPM = EWPM IF (NSPM > MaxPM) MaxPM = NSPM IF (EWPM < MinPM) MinPM = EWPM IF (NSPM < MinPM) MinPM = NSPM IF (EWPMErr > MaxPMErr) MaxPMErr = EWPMErr IF (NSPMErr > MaxPMErr) MaxPMErr = NSPMErr IF (EWPMErr < MinPMErr) MinPMErr = EWPMErr IF (NSPMErr < MinPMErr) MinPMErr = NSPMErr IF (EpochRA > MaxEpoch) MaxEpoch = EpochRA IF (EpochDec > MaxEpoch) MaxEpoch = EpochDec IF (EpochRA < MinEpoch) MinEpoch = EpochRA IF (EpochDec < MinEpoch) MinEpoch = EpochDec !----------------------------------------------------------------------------------------------------------------------------------! ! Decide what to do with the magnitude information ! ! ! ! B magnitude histogram shows distribution peaks at 20.2, drops to half the peak at 21.0, with a big drop after 22.0 ! ! R magnitude histogram shows distribution peaks at 18.7, drops to half the peak at 19.9, with a big drop after 21.0 ! ! ! ! Don't trust any B fainter than 22 or brighter than -2 (assign a large default) ! ! Don't trust any R fainter than 21 or brighter than -2 (assign a large default) ! ! ! ! Consider using the JHK magnitudes, where available, to derive color information that could be used to flag bad USNO magnitudes ! !----------------------------------------------------------------------------------------------------------------------------------! IF (B1Mag > 22.0 .OR. B1Mag < -2.0) B1Mag = 99.99 IF (B2Mag > 22.0 .OR. B2Mag < -2.0) B2Mag = 99.99 IF (R1Mag > 21.0 .OR. R1Mag < -2.0) R1Mag = 99.99 IF (R2Mag > 21.0 .OR. R2Mag < -2.0) R2Mag = 99.99 ! If a source has no believable first or second epoch B or R magnitude, increment the "no magnitude information" counter IF (B1Mag > 50.0 .AND. B2Mag > 50.0 .AND. R1Mag > 50.0 .AND. R2Mag > 50.0) NumNoMags = NumNoMags + 1 ! If one of the two magnitudes changed by more than a mangitude between epochs, flag the source as variable and increment the ! variable counter; don't treat as variable those cases where one of the magnitudes has been assigned a large default for either ! not having any magnitude information of for not being believable by ignoring differences larger than 100 BDiff = ABS(B1Mag - B2Mag) RDiff = ABS(R1Mag - R2Mag) IF ((BDiff > 1.0 .AND. BDiff < 50.0) .OR. (RDiff > 1.0 .AND. RDiff < 50.0)) THEN Variable = .TRUE. NumVar = NumVar + 1 ELSE Variable = .FALSE. END IF ! If the magnitude difference is too large, pick the brighter of the two; note that this is an arbitrary choice and could ! be revisited; if the source isn't wildly variable, use the average of the magnitudes at the two epochs ! Note that a threshold of 1 mag represents very roughly 3 sigma, given that the photographic photometry is typically good ! to only 0.3 mag IF (BDiff > 1.0) THEN BMag = MIN(B1Mag,B2Mag) ELSE BMag = (B1Mag + B2Mag)/2.0 END IF IF (RDiff > 1.0) THEN RMag = MIN(R1Mag,R2Mag) ELSE RMag = (R1Mag + R2Mag)/2.0 END IF !----------------------------------------------------------------------------------------------------------------------------------! ! Section for flagging anomalies and saving them to a file ! !----------------------------------------------------------------------------------------------------------------------------------! ! Flag stars with magnitudes between 0 and -2 (should be very few) IF (B1Mag > -2.0 .AND. B1Mag < 0.0) WRITE (198,'(2A,F6.2)') Line(1:35),' B1 ',B1Mag IF (B2Mag > -2.0 .AND. B2Mag < 0.0) WRITE (198,'(2A,F6.2)') Line(1:35),' B2 ',B2Mag IF (R1Mag > -2.0 .AND. R1Mag < 0.0) WRITE (198,'(2A,F6.2)') Line(1:35),' R1 ',R1Mag IF (R2Mag > -2.0 .AND. R2Mag < 0.0) WRITE (198,'(2A,F6.2)') Line(1:35),' R2 ',R2Mag ! Flag stars with proper motions in excess of 5 arcsec per year (should be very few) IF (ABS(EWPM*3600) > 5.0) THEN NumHighPM = NumHighPM + 1 WRITE (198,'(2A,F6.2,A)') Line(1:35),' east-west proper motion ',EWPM*3600,' arcsec/yr' END IF IF (ABS(NSPM*3600) > 5.0) THEN NumHighPM = NumHighPM + 1 WRITE (198,'(2A,F6.2,A)') Line(1:35),' north-south proper motion ',NSPM*3600,' arcsec/yr' END IF ! Flag stars with mean epoch differences IF (ABS(EpochRA - EpochDec) > 0.005) THEN NumEpochD = NumEpochD + 1 WRITE (198,'(2A,2F8.2)') Line(1:35),' epoch difference',EpochRA,EpochDec END IF !----------------------------------------------------------------------------------------------------------------------------------! ! Encode the various catalog values to be saved ! ! ! ! Right ascensions and declinations are given to millionths of a degree precision, which is 3.6 milliarcseconds, so 1 milliarcsec ! ! precision is enough to prevent loss of precision while not overflowing a 32-bit integer value ! ! ! ! Right ascensions will range from 0 to 1,296,000,000 milliarcsec, well under the 2**31 - 1 limit for a four-byte integer ! ! Declinations will range from -324,000,000 to +324,000,000, well under the 2**31 - 1 limit for a four-byte integer ! ! ! ! Position errors range from 1 to 210 milliarcsec, thus can fit in a single byte with milliarcsec precision if biased by -100 ! ! ! ! Proper motions errors are as small as 0.1 milliarcsec per year, and proper motions are as big as 10328 milliarcsec per year, so ! ! they won't fit in a two-byte integer without loss of precision; need to use four-byte integers, even though overkill ! ! ! ! Proper motion errors are as big as 89.3 milliarcsec per year; need to use two-byte integers, even though overkill ! ! ! ! Epochs are given to hundredths of a year and span 1909 to 2006, so subtract 1900 from all both RA and DEC epochs, then multiply ! ! by 100; encoded values should range from 969 to 10516, which can be handled in a two-byte integer ! ! ! ! B and R magnitudes are given to 0.01 mag precision (mostly unjustified); could be handled in a single byte if precision were ! ! reduced to 0.1 mag, but otherwise needs two-byte integers, as long as the large default is less than 327.67 mag ! !----------------------------------------------------------------------------------------------------------------------------------! IRA = NINT(RA *DegToMAS) IDec = NINT(Dec*DegToMAS) IEWErr = NINT(EWErr*DegToMAS) - 100 INSErr = NINT(NSErr*DegToMAS) - 100 IEWPM = NINT(EWPM*DegToMAS*10) INSPM = NINT(NSPM*DegToMAS*10) IEWPMErr = NINT(EWPMErr*DegToMAS*10) INSPMErr = NINT(NSPMErr*DegToMAS*10) IEpochDec = NINT((EpochDec - 1900.0)*100) IEpochRA = NINT((EpochRA - 1900.0)*100) IBMag = NINT(BMag*100) IRMag = NINT(RMag*100) ! Write record to the appropriate declination zone and increment the zone population counter WRITE (ZoneNum,'(I10,I11,2I4,2I7,2I4,2I6,2I5)') IRA,IDec,IEWErr,INSErr,IEWPM,INSPM,IEWPMErr,INSPMErr,IEpochRA,IEpochDec,IBMag, & IRMag ZoneCount(ZoneNum) = ZoneCount(ZoneNum) + 1 ! Progress indicator; note that most systems default to logical unit number 6 for standard output, but in this program, that ! unit number is being used by a declination zone, so writing to the standard output unit may fail, unless the compiler is smart ! enough to preconnect a different unit number; otherwise, explicitly open an unused unit like 200 using the CON device (Windows) ! or /dev/stdout (Unix) and replace the asterisks with that unit Progress = Progress + 1 IF (MOD(Progress,100000) .EQ. 0) WRITE (*,*) Progress END DO WRITE (*,'(I9,A)') Progress,' finished reading catalog' !----------------------------------------------------------------------------------------------------------------------------------! ! Okay, reading of the monster catalog is now complete ! ! ! ! Save some statistics to help decide on an encoding algorithm ! !----------------------------------------------------------------------------------------------------------------------------------! WRITE (197,'(A,I0)') 'Total number of records = ',Progress WRITE (197,'(A,F9.6,A)') 'Maximum error = ',3600*MaxErr,' arcsec' WRITE (197,'(A,F9.6,A)') 'Minimum error = ',3600*MinErr,' arcsec' WRITE (197,'(A,F9.6,A)') 'Maximum PM = ',3600*MaxPM,' arcsec/yr' WRITE (197,'(A,F9.6,A)') 'Minimum PM = ',3600*MinPM,' arcsec/yr' WRITE (197,'(A,F9.6,A)') 'Maximum PM error = ',3600*MaxPMErr,' arcsec/yr' WRITE (197,'(A,F9.6,A)') 'Minimum PM error = ',3600*MinPMErr,' arcsec/yr' WRITE (197,'(A,F7.2)') 'Maximum mean epoch = ',MaxEpoch WRITE (197,'(A,F7.2)') 'Minimum mean epoch = ',MinEpoch WRITE (197,'(A,I0)') 'Number of magnitudes = ',NumMags WRITE (197,'(A,I0)') 'Number with no mags = ',NumNoMags WRITE (197,'(A,F5.2)') 'Faintest source = ',MaxMag WRITE (197,'(A,F6.2)') 'Brightest source = ',MinMag WRITE (197,'(A,I0)') 'Number variables = ',NumVar WRITE (197,'(A,I0)') 'Number high prop motion = ',NumHighPM WRITE (197,'(A,I0)') 'Number epoch difference = ',NumEpochD WRITE (197,'(A)') 'Number of sources per declination zone' ! Save number of entries per declination zone, and while we're at it, close each declination zone file, freeing those logical ! unit numbers Sum = 0 DO I=0,179,1 WRITE (197,'(I3.3,I9)') I,ZoneCount(I) CLOSE (I) Sum = Sum + ZoneCount(I) END DO WRITE (197,'(A,I10,A)') 'Sum of zones counts = ',Sum,' (should be same as total number of records)' ! Save magnitude histograms DO I=-110,670,1 WRITE (196,'(I4,2I11)') I,BMagArray(I),RMagArray(I) END DO !----------------------------------------------------------------------------------------------------------------------------------! ! PHASE 2 - Sort each declination zone by right ascension ! ! ! ! Note that this section also writes to standard output; see the note above just before the progress indicator increment ! !----------------------------------------------------------------------------------------------------------------------------------! DO I=0,179,1 ! Open raw, sorted catalog, and accelerator files ! Note that FORM='BINARY' is a Compaq Visual Fortran extension; change to ACCESS='STREAM' for a compiler that supports the ! Fortran 2003/2008 standards WRITE (File,'(I3.3,A4)') I,'.raw' OPEN (1,FILE=File) WRITE (File,'(I3.3,A4)') I,'.cat' OPEN (2,FILE=File,FORM='BINARY') WRITE (File,'(I3.3,A4)') I,'.acc' OPEN (3,FILE=File) ! Allocate memory for arrays Size = ZoneCount(I) ALLOCATE (RAArray(Size)) ALLOCATE (Ptr(Size)) ALLOCATE (Rest(Size)) ! Read in all entries in a declination zone DO J=1,Size,1 READ (1,'(I10,A63)',IOSTAT=ErrorCode) RAArray(J),Rest(J) IF (ErrorCode < 0) THEN WRITE (*,'(A,I3.3,A)') 'Premature end of file encountered for ',I,'.raw' WRITE (197,'(A,I3.3,A)') 'Premature end of file encountered for ',I,'.raw' EXIT END IF END DO ! Sort by right ascension WRITE (*,'(A,I3.3)') 'Starting sort of zone ',I CALL SortIntegerShell(RAArray,Size,Ptr) WRITE (*,'(A)') 'Finished' ! Write sorted array and generate accelerator table with an entry every degree; append zone size to accelerator file ! to facilitate end-of-file detection when using direct access methods (which should only trigger an error condition, ! not an end-of-file condition) Limit = 0 DO J=1,Size,1 READ (Rest(Ptr(J)),'(I11,2I4,2I7,2I4,2I6,2I5)') IDec,IEWErr,INSErr,IEWPM,INSPM,IEWPMErr,INSPMErr,IEpochRA,IEpochDec,IBMag,IRMag WRITE (2) RAArray(Ptr(J)),IDec,IEWErr,INSErr,IEWPM,INSPM,IEWPMErr,INSPMErr,IEpochRA,IEpochDec,IBMag,IRMag IF (RAArray(Ptr(J)) > Limit) THEN WRITE (3,'(I8)') J Limit = Limit + DegToMAS END IF END DO WRITE (3,'(I8)') Size ! Free memory DEALLOCATE (RAArray) DEALLOCATE (Ptr) DEALLOCATE (Rest) ! Close files CLOSE (1) CLOSE (2) CLOSE (3) END DO !----------------------------------------------------------------------------------------------------------------------------------! ! SUBROUTINE SortIntegerShell ! ! Sorts an integer array indirectly via a pointer array (Shell-Mezgar algorithm) ! ! ! ! Copyright (C) 1999 by David J. Tholen ! ! Last update: 1999 Mar 29 ! !----------------------------------------------------------------------------------------------------------------------------------! CONTAINS SUBROUTINE SortIntegerShell(IA,N,PT) IMPLICIT NONE INTEGER I INTEGER IA (:) INTEGER J INTEGER K INTEGER L INTEGER M INTEGER N INTEGER PT (:) INTEGER TMP !----------------------------------------------------------------------------------------------------------------------------------! DO I=1,N,1 PT(I) = I ! initialize the pointer array END DO ! The variable M contains the distance between the two array elements to be compared; the algorithm starts with this distance ! as half the array size, and becomes progressively smaller by a factor of two; the sort is completed after adjacent elements ! (M = 1) have been compared M = N DO M = M/2 IF (M .LE. 0) RETURN L = N - M ! work up through array DO K=1,L,1 I = K ! Initialize indicies for comparison 3 CONTINUE J = I + M IF (IA(PT(I)) .GT. IA(PT(J))) THEN TMP = PT(I) PT(I) = PT(J) PT(J) = TMP I = I - M IF (I .GE. 1) GO TO 3 END IF END DO END DO ! loop back to divide distance in half again END SUBROUTINE SortIntegerShell END PROGRAM PPMXL