FTN77 $FILES 0,1 PROGRAM DSCN8 (3,70) ****************************************************************** * in file /MIKE/PROGS/DSCN8.FTN on the control computer * version <901109.1011> * For drift scans through extended sources at any frequency. * Uses disc source file for input data * Run string used for LU, file name, start source number * * Programs called are : * DPRC2 precess 1950 RA,DEC to current epoch * DRIVE produce drive command for COMND * DSCPL plot drift scan * ALARM set servo alarm on if fault occurs * files INCLUDED - see end of this listing * PARAMETER (MAXI = 6, MAXDATA = 480, NSPBW = 12) LOGICAL FILEND, START, FAULT, GFIT, IFBRK, & PLOTREPEAT, DRIVEWAIT, TRACK, NOTEND CHARACTER NAME*20, FILENAME*10, CLU*2, CNSORCE*3, & SYSTEM(MAXI)*4 INTEGER*2 IANG(5), ITIME(5), IPRAM(5), ISTORE(1006), ANTBUF(7), & IDVMSL(MAXI), NDSEL(MAXI), ISYNTH(MAXI), DRIVE(3) REAL DATA(MAXDATA), AVRG(MAXDATA), & RABMOFF(MAXI), DCBMOFF(MAXI), DNOISE(6), WAVE(MAXI), & HPBW(MAXI), BWFN(MAXI), RADIOMETER(MAXI) REAL*8 DFREQ * EQUIVALENCE * data buffer for DPRC2 & (IANG(1),RAP), (IANG(3),DECP), * data buffer for DSCPL & (ISTORE(1),DATA(1)), (ISTORE(961),NAME(1:1)), & (ISTORE(971),WAVE(1)), (ISTORE(983),IYR), & (ISTORE(984),ITIME(1)), (ISTORE(989),DEC), & (ISTORE(991),RAS50), (ISTORE(993),RAE50), & (ISTORE(995),LAST), (ISTORE(996),IEPOCH), & (ISTORE(997),INDEX), (ISTORE(998),LUOUT), & (ISTORE(999),LP), (ISTORE(1000),IXPOS), & (ISTORE(1001),IYPOS), (ISTORE(1002),IXC), & (ISTORE(1003),IYC), (ISTORE(1004),IHITE), & (ISTORE(1005),IWIDE), (ISTORE(1006),MKNUM), * data buffer for DRIVE & (ANTBUF(1), ITYPE), (ANTBUF(2), RADRIVE), & (ANTBUF(4), DECDRIVE), (ANTBUF(6), TRACK), & (ANTBUF(7), DRIVEWAIT) * COMMON /LUS/ LU,MTLU,LU3330,LU3336,LU8662,LU8840A,LU8840B * * index 1 2 3 4 5 6 DATA SYSTEM /'2.5S', '3.6D', ' 6S', ' 6D', ' 13S', ' 18S' /, & WAVE / 2.5, 3.6, 6.1, 6.2, 13.0, 18.0 /, & HPBW / 0.066, 0.096, 0.160, 0.160, 0.333, 0.493 /, * exptl hpbw's 0.106 0.163 0.338 1990 jun & BWFN / 0.150, 0.450, 0.470, 0.830, 0.620, 1.250 /, & RABMOFF /-0.405, +0.154, +0.555, +0.064, +0.057, +0.056 /, & DCBMOFF /-0.050, -0.565, -0.162, -0.181, +0.445, +1.700 /, & DNOISE / 148.0, 10.0, 14.0, 1.54, 3.96, 11.2 /, & IDVMSL / 77B, 11b, 33b, 33b, 55b, 77b /, & NDSEL / 2, 1, 2, 3, 4, 5 /, & RADIOMETER/ 18.0, 3.6, 6.0, 6.0, 13.0, 18.0 /, & ISYNTH / 8662, 3330, 8660, 8660, 5100, 8662 /, * & LP /12/, LUDVM/35/, IFF/14b/, FILENAME /' '/, & DTUS /239.33/, DTSS /240.0/, DTR /0.017453293/, & NTSETTLE /-3/, & DRIVE / 'DRIVE '/, IPAGE / 1 / * * WAVE nominal operating wavelength (cm) * HPBW half power beam width (degrees) * BWFN beam width to first null (degrees) * RABMOFF beam offset in right ascension (degrees) * DCBMOFF beam offset from axis in declination (degrees) * DNOISE noise tube calibration (Kelvins) * IDVMSL dvm select code for unfiltered output * NDSEL calibration noise diode selection code * RADIOMETER radiometer selection * LU,LP I/O and Plotter logical units * LUDVM lu of dvm in use (35 or 36) * DTUS degrees RA to universal seconds conversion * DTSS degrees RA to sidereal seconds " * DTR degrees to radians " * MAXDATA data array size for drift scan * NSPBW number of data points per half power beamwidth * NTSETTLE noise cal settling time (seconds) * IPAGE single page chart advance on plotter * * Initialise the program * ---------------------- * * pickup i/o lu, filename and start source number from run string CALL FPARM (CLU,FILENAME,CNSORCE) * * set up user terminal LU READ (CLU,*,ERR=20) LU 20 IF (LU .NE. 1 .AND. LU .NE. 7) LU = 7 * set standard input, output device for HP 1000 CALL FSYSU (LU,LU) LUOUT = LU * do a form feed on the Anadex printer * IF (LU .EQ. 7) PRINT *, CHAR(14b) * READ (CNSORCE,*,ERR=30) NSORCE 30 IF (NSORCE .LT. 1) NSORCE = 1 * * clear the switch register on the HP 1000 CALL ISSR(0) START = .TRUE. FAULT = .FALSE. * CALL EXEC(11,ITIME,IYR) WRITE (LU,1000) IYR, ITIME(5) 1000 FORMAT( 3X,'DSCN8',I5,' DAY',I3/3X,17'='// * 3X,'Set S reg. bit 13 on to suspend'/ * 18X,'14 on for next source'/ * 18X,'15 on to end.'/) * * start point for each source * --------------------------- * * Read the source parameters from the disc file 100 CALL DSCRD (LU,START,SYSTEM,FILENAME,NSORCE,FILEND, * NAME,IEPOCH,RA,DEC,SCANL,INDEX, * NSCANS,GFIT,PLOTREPEAT,DFREQ) IF (FILEND) THEN GO TO 900 END IF START = .FALSE. * * pick up latest beam offsets from SYPARM::16 CALL RDSYS (LU,INDEX,SYSTEM,RABMOFF,DCBMOFF) * * Set up parameters for the drift scan * ------------------------------------ * * default to precessed coordinates RAP = RA DECP = DEC IF (IEPOCH .EQ. 1950) THEN * precess source coordinates CALL EXEC(9,5HDPRC2,IANG(1),IANG(2),IANG(3),IANG(4)) CALL RMPAR(IANG) END IF * * cosine of precessed declination COSDECP = COS (DECP * DTR) * Calculate beamwidths in degrees of RA at source declination HPBWDEC = HPBW(INDEX) / COSDECP BWFNDEC = BWFN(INDEX) / COSDECP * * maximum scan length SCANMAX = MAXDATA / NSPBW * HPBWDEC * compare with wanted scan length SCANL = MIN (SCANL, SCANMAX) * minimum scan length IF (GFIT) THEN * must be long enough to fit baseline for Gaussian fit SCANMIN = BWFNDEC + 2 * HPBWDEC ELSE * must be long enough to include first minima of beam SCANMIN = BWFNDEC END IF SCANL = MAX (SCANL, SCANMIN) * * period for NSAM samples (seconds) SAMINT = HPBWDEC * DTUS / NSPBW * * Find current RA and DEC with beam offsets added ONRA = RAP + RABMOFF(INDEX) / COSDECP ONDEC = DECP + DCBMOFF(INDEX) * * Find start and end RA in current and 1950 coordinates HALFSC = SCANL / 2.0 RASTRT = POSRA (RAP - HALFSC) RAEND = POSRA (RAP + HALFSC) RAS50 = POSRA (RA - HALFSC) RAE50 = POSRA (RA + HALFSC) DELTRA = RAEND - RAE50 * * plot number on current page & number of scans averaged NPLOT = 0 ADAV = 0.0 * PRINT *, CHAR(14B) WRITE (LU,1700) & NSORCE, NAME, & RA, DEC, IEPOCH, SCANL, SYSTEM(INDEX), NSCANS, IPAGE, & RABMOFF(INDEX), DCBMOFF(INDEX), RAP, DECP, DNOISE(INDEX), & SAMINT 1700 FORMAT & (/3X,'SOURCE',I4,' : ',A, & /6X,'RA',6X,'DEC',4X,'epoch',2X,'length',2X,'system', & 2X,'scans',2X,'pages' & /3X,2F8.3,I6,F6.2,'deg',A7,4X,i3,4X,I2, & //6X,'beam offsets',7X,'RAP',5X,'DECP',2X,'noise cal.', & 2X,'sample period', & /3X,2F8.3,3X,2F8.3,F8.2,'K',6X,F6.2,'S') * IF (DFREQ .GT. 0D0) THEN PRINT *,' Set LO synthesizer',ISYNTH(INDEX), & ' to',DFREQ,'MHz' CALL SETLO (ISYNTH,INDEX,DFREQ) ELSE PRINT *,' LO synthesizer not set' END IF * * Start of loop for each drift scan * --------------------------------- * DO LOOP = 1,NSCANS WRITE (LU,'(/3x,''SCAN '',i3,'' : _'')') LOOP * * reset the "computer failed" alarm CALL OBSOK * * plot control NPLOT = NPLOT + 1 IF (NPLOT .GT. 4) NPLOT = 1 * select DVM inputs INSEL = IDVMSL(INDEX) CALL EXEC (2,113B,INSEL,1,5,5) * * Find command HA, compare to east horizon limit 410 COMHA = TIME(2) / DTSS - RASTRT COMHA = PMOD(COMHA) WRITE (LU,4100) COMHA 4100 FORMAT(3X,'Command HA = ',F7.3) * * West limit is west horizon * East limit is one hour below east horizon IF (COMHA .GT. HALIM(+90.0,ONDEC) .OR. * COMHA .LT. (HALIM(-90.0,ONDEC) -15.0)) THEN * Source will not rise in less than 1 hour - try next source WRITE (LU,4120) 4120 FORMAT(6X,'Source > 1 hour below horizon') GO TO 800 END IF * * Else calculate wait time until source rises, -ve for wait 415 ISWAIT = (COMHA - HALIM(-90.0,ONDEC)) * DTUS IF (ISWAIT .LT. -1 .AND. COMHA .GT. 0.0) THEN * Source below eastern horizon, wait until it has risen MWAIT = -ISWAIT / 60 WRITE (LU,4150) MWAIT 4150 FORMAT(3X,'WAIT',I5,' M') CALL EXEC(12,0,2,0,ISWAIT) END IF * * Drive to start coordinates * -------------------------- * allow for calibration time, add beam offsets 420 RADRIVE = POSRA (RASTRT - 0.30 * + RABMOFF(INDEX) / COSDECP) DECDRIVE = DECP + DCBMOFF(INDEX) ITYPE = 4 TRACK = .FALSE. DRIVEWAIT = .TRUE. * call drive program to send command to STEER CALL EXEC (9,DRIVE,0,0,0,0,0,ANTBUF,7) CALL EXEC (9,DRIVE,0,0,0,0,0,ANTBUF,7) * * Do a noise tube calibration, return dvm units/K conversion CALL CALT2 (LU,LUDVM,NDSEL,INDEX,NTSETTLE,DNOISE, & GASCAL,FAULT) IF (FAULT) THEN * abort if repeated calibration failures GO TO 900 END IF * * find current antenna position CALL ANTPOS (RELHA,RELDEC) * * Western horizon check, end scan if failed 425 IF (RELHA .GT. (HALIM(+90.0,RELDEC)-0.5)) GO TO 700 * * Wait for beginning of drift scan * * antenna HA + Beam offset + Pointing correction TRUHA = RELHA + RABMOFF(INDEX)/COS(RELDEC*DTR) & + HAPC(RELHA,RELDEC) * true current RA of telescope RANOW = RATEL(TRUHA) WRITE (LU,4200) RELHA, RELDEC, RAS50, RAE50, RASTRT, & RAEND, RANOW 4200 FORMAT(/4X,'RELHA',5X,'RELDEC',4X,'RAS50',5X,'RAE50', & 5X,'RASTRT',4X,'RAEND',5X,'RANOW'/ & 7(F10.3)) * RAENDMOD = RAEND * 360 - 0 degree RA continuity IF (RASTRT .GT. 180.0 .AND. RAEND .LT. 180.0) & RAENDMOD = RAEND + 360.0 * * time to wait in seconds ISWAIT = -((RASTRT - RANOW) * DTUS - SAMINT / 2) WRITE (LU,4300) -ISWAIT 4300 FORMAT(3X,'wait',I5,' S_') * * If start time has passed, start again IF (ISWAIT .GT. -1) THEN GO TO 420 ELSE CALL EXEC (12,0,2,0,ISWAIT) END IF * * Start of drift scan * ------------------- * CALL EXEC (11,ITIME,IYEAR) WRITE (LU,'(6x,''Start scan at '',4i5,'' _'')') & (itime(i) , i = 5,2,-1) * DO I = 1,MAXDATA * initialize data point, data counter, value sum DATA(I) = 0 COUNTS = 0 END = TIME(1) + SAMINT NOTEND = .TRUE. * DO WHILE (NOTEND) * read from fluke multimeter READ (LUDVM,*) VALUE DATA(I) = DATA(I) + VALUE COUNTS = COUNTS + 1 NOTEND = (TIME(1) .LT. END) END DO * * take mean value DATA(I) = DATA(I) / COUNTS * * increment last point counter LAST = I * * end data capture if end RA reached RANOW = RATEL(TRUHA) IF (RASTRT .GT. 180.0 .AND. RAEND .LT. 180.0 .AND. * RANOW .LT. 180.0) RANOW = RANOW + 360.0 IF (RANOW .GE. RAENDMOD) go to 600 END DO * * If data array full, find end RA 600 RAEND = RATEL(TRUHA) - SAMINT / (DTSS * 2) RAE50 = POSRA (RAEND - DELTRA) * * end of drift scan - scale data to Kelvins, add to average * --------------------------------------------------------- WRITE (LU,'(6X,''End scan''/)') * DO I = 1,MAXDATA IF (I .GT. LAST) DATA(I) = DATA(LAST) DATA(I) = DATA(I) / GASCAL AVRG(I) = (AVRG(I)*ADAV + DATA(I)) / (ADAV + 1.0) END DO * increment number of averaged spectra ADAV = ADAV + 1.0 IADAV = ADAV * * do a gaussian fit if wanted IF (GFIT) THEN WRITE (LU,'(6x,''this scan :'')') CALL DOGFT (LU,LAST,DATA,SYSTEM,INDEX) IF (IADAV .GE. 2) THEN WRITE (LU,'(6x,''average scan :'')') CALL DOGFT (LU,LAST,AVRG,SYSTEM,INDEX) END IF END IF * * Check switch registers and number of scans for plot options * ----------------------------------------------------------- * * Plot individual scans if number done is greater than 1, * or bits 14 and 15 are clear and more than 1 scan wanted IF (IADAV .GT. 1 .OR. (ISSW(14) .GE. 0 .AND. * ISSW(15) .GE. 0 .AND. * NSCANS .gt. 1)) THEN * set position and character size for individual plots IXPOS = 1000 IYPOS = 2125 + (NPLOT - 1) * 2250 IXC = 100 IYC = 150 IHITE = 2100 IWIDE = 4000 - IYC * 3 / 2 MKNUM = 3 * * Call DSCPL to plot individual scan IPL = 1 IF (NPLOT .GT. 1) IPL = 2 C ICLAS = 0 C CALL CLRQ(1,ICLAS) C CALL EXEC(20+140000B,0,ISTORE,1006,0,0,ICLAS,*890) C CALL EXEC(24,5HDSCPL,ICLAS,IPL,IPAGE) CALL EXEC(24,5HDSCPL,IPL,IPAGE,0,0,0,ISTORE,1006) * * plot average scan DO IA = 1,MAXDATA DATA(IA) = AVRG(IA) END DO IXPOS = 5990 C ICLAS = 0 C CALL CLRQ(1,ICLAS) C CALL EXEC(20+140000B,0,ISTORE,1006,0,0,ICLAS,*890) C CALL EXEC(24,5HDSCPL,ICLAS,IPL,0) CALL EXEC(24,5HDSCPL,IPL,0,0,0,0,ISTORE,1006) END IF * * If next source, or end bits set, or last scan done IF (ISSW(14) .LT. 0 .OR. ISSW(15) .LT. 0 * .OR. IADAV .EQ. NSCANS) THEN IF (IADAV .EQ. 1 * .OR. (IADAV .GT. 1 .AND. PLOTREPEAT)) THEN * plot full page average scan DO IA = 1,MAXDATA DATA(IA) = AVRG(IA) END DO IXPOS = 1500 IYPOS = 5500 IXC = 150 IYC = 180 IHITE = 8000 IWIDE = 8500 - IYC * 3 / 2 MKNUM = 7 C ICLAS = 0 C CALL CLRQ(1,ICLAS) C CALL EXEC(20+140000B,0,ISTORE,1006,0,0,ICLAS,*890) C CALL EXEC(24,5HDSCPL,ICLAS,3,IPAGE) CALL EXEC(24,5HDSCPL,3,IPAGE,0,0,0,ISTORE,1006) END IF * Print average scan IF ( .NOT. GFIT) THEN CALL DNORM (DATA,LAST,WAVE,INDEX) END IF CALL TABLE (LU,LAST,SCANL,DATA) endif * * Check for program suspension IF (ISSW(13) .LT. 0 .OR. IFBRK(0)) THEN PAUSE 'Suspended - type *GO,DSCN8 to restart' CALL ISSR (0) WRITE (LU,'(''Set S-reg. options now for DSCN8'')') CALL EXEC (12,0,2,0,-30) END IF * * Check for end IF (ISSW(15) .LT. 0 .OR. FAULT) GO TO 900 * * Check for next source IF (ISSW(14) .LT. 0) GO TO 800 * 700 END DO * * Go to next source * ----------------- * 800 NSORCE = NSORCE + 1 GO TO 100 * * fault calling DSCPL C 890 WRITE (LU,'(''error return from class write to DSCPL'')') C FAULT = .TRUE. * * end * --- * 900 IF (FAULT) CALL EXEC(10,5HALARM) STOP 'Bye from DSCN8' END * * SUBROUTINE DSCRD (LU,START,SYSTEM,FILENAME,NSORCE, FILEND, & NAME,IEPOCH,RA,DEC,SCANL,INDEX,NSCANS,GFIT,PLOTREPEAT,DFREQ) ************************************************************************ * Read source data from disc file * inputs : * LU integer*2 i/o lu * START logical*2 at start of program or not * SYSTEM character*4 array of system designations * FILENAME character*10 source file name * NSORCE integer*2 number of source in file to read * outputs : * FILEND logical*2 true if end of source file encountered * NAME character*20 name identifying the source * IEPOCH integer*2 year of coordinates (e.g. 1950, 1990) * RA real*4 RA in decimal degrees * DEC real*4 DEC in decimal degrees * SCANL real*4 drift scan length in degrees of RA * INDEX integer*2 internal array index identifying the system * NSCANS integer*2 number of drift scans to do on this source * GFIT logical*2 true for Gaussian fit, else false * PLOTREPEAT logical*2 true for repeat of final average plot * DFREQ real*8 LO syynthesizer frequency too be set * File format : * NAME 20 character name, in quotes * SYSID 4 character system identifier, in quotes * IEPOCH integer coordinate epoch - 1950 or current year * RA real Right Ascension in decimal degrees * DEC real Declination in decimal degrees * SCANL real Wanted scan length in degrees of RA * NSCANS integer Number of drift scans to be done on the source * IGFIT integer 1 to do a gaussian fit to source, else 0 * IREPEAT integer 1 for repeat of final average plot, else 0 * DFREQ real*8 L O synthesizer frequency to be set * PARAMETER (MAXI = 6) CHARACTER FILENAME*10, NAME*20, CBUFFER*80, SYSTEM(MAXI)*4, & SYSID*4, SYSTEMP*4 LOGICAL START, FILEND, DATAERROR, GFIT, PLOTREPEAT REAL*8 DFREQ * * Request file info if starting * ----------------------------- IF (START .AND. (FILENAME(1:1) .EQ. ' ' .OR. NSORCE .LE. 0)) THEN 110 WRITE (lu,'(3X,''Disc file name ? _'')') READ (lu,'(a)',err=110) filename(1:10) 120 WRITE (lu,'(/3X,''Start at source number ? _'')') READ (lu,*,err=120) nsorce END IF * * defaults : NSCANS = 1 IGFIT = 0 IREPEAT = 0 FILEND = .FALSE. DFREQ = 0D0 * OPEN (99,FILE=FILENAME,IOSTAT=IOS,ERR=700) * 200 REWIND (99) DO 210 I = 1,NSORCE READ (99,'(A)',IOSTAT=IOS,ERR=800,END=900) CBUFFER 210 END DO * * truncate read at last parameter in list in character buffer CBUFFER(79:79) = '/' READ (CBUFFER,*,IOSTAT=IOS,ERR=800,END=300) & NAME,SYSID,IEPOCH,RA,DEC,SCANL,NSCANS,IGFIT,IREPEAT,DFREQ * CLOSE (99) * * Check correctable data 300 IF (DEC .GT. 270.0) DEC = DEC - 360.0 IF (NSCANS .LT. 1) NSCANS = 1 GFIT = .FALSE. IF (IGFIT .EQ. 1) GFIT = .TRUE. PLOTREPEAT = .FALSE. IF (IREPEAT .EQ. 1) PLOTREPEAT = .TRUE. CALL UPCASE (NAME) CALL UPCASE (SYSID) * * Check for fatal data errors DATAERROR = .FALSE. * * right justify the input system index DO WHILE (SYSID(4:4) .EQ. ' ') SYSTEMP = ' ' SYSTEMP(2:4) = SYSID(1:3) SYSID = SYSTEMP END DO * * find the system index INDEX = 0 DO I = 1, MAXI IF (SYSTEM(I) .EQ. SYSID) INDEX = I END DO IF (INDEX .EQ. 0) THEN DATAERROR = .TRUE. END IF * 320 IF (IEPOCH .NE. 1950 .AND. IEPOCH .LT. 1990 * .OR. RA .LT. 0.0 .OR. RA .GT. 360.0 * .OR. DEC .LT. -90.0 .OR. DEC .GT. 45.0 * .OR. SCANL .LE. 0.0) * DATAERROR = .TRUE. * IF ( .NOT. DATAERROR) THEN RETURN ELSE * Fatal error in data GO TO 800 END IF * * error opening file 700 WRITE (lu,'(''Error '',i4,'' opening file '',a10)') * ios, filename STOP 'BYE' * * read error encountered, try next source 800 WRITE (lu,'(''Error'',i5,'' at source'',i3/a)') & ios, nsorce, cbuffer(1:80) NSORCE = NSORCE + 1 GO TO 200 * * end of file 900 WRITE (lu,'(''End of file at source'',i4)') nsorce FILEND = .TRUE. RETURN END * SUBROUTINE DOGFT (LU,LAST,ARRAY,SYSTEM,INDEX) ************************************************************************ * set up Gaussian fit simulating the beam * inputs : * LU integer*2 i/o lu * LAST integer*2 last valid point in the array * ARRAY real*4 data array * SYSTEM character*4 array of system identifiers * INDEX integer*2 index identifying system in use * PARAMETER (MAXI = 6) LOGICAL*2 DUALBEAM REAL*4 ARRAY(1), FLUXESTA(11), FLUXESTB(11) CHARACTER SYSTEM(MAXI)*4 * * subtract linear fit to baseline of array CALL LNOFF (LAST,ARRAY,RMS) WRITE (lu,'(9x,''rms noise = '',1pg12.6,'' K'')') RMS * * define the mid beam points in the data array IF (SYSTEM(INDEX) .EQ. ' 6D') THEN MIDBEAMA = LAST / 2 - 8 MIDBEAMB = LAST / 2 + 8 DUALBEAM = .TRUE. ELSE IF (SYSTEM(INDEX) .EQ. '3.6D') THEN MIDBEAMA = LAST / 2 - 13 MIDBEAMB = LAST / 2 + 13 DUALBEAM = .TRUE. ELSE * for single feeds MIDBEAM = LAST / 2 DUALBEAM = .FALSE. END IF * IF (DUALBEAM) THEN * make 11 estimates of flux for each beam * error in flux is rms/2 per beam owing to gaussian fit J = 0 DO I = MIDBEAMA - 5, MIDBEAMA + 5 J = J + 1 CALL GFLUX (ARRAY,I,FLUXESTA(J)) END DO J = 0 DO I = MIDBEAMB - 5, MIDBEAMB + 5 J = J + 1 CALL GFLUX (ARRAY,I,FLUXESTB(J)) END DO * * find (signed) maximum flux in each beam FLUXA = FLUXESTA(1) FLUXB = FLUXESTB(1) DO J = 2, 11 FLUXA = MAX (FLUXA, FLUXESTA(J)) FLUXB = MIN (FLUXB, FLUXESTB(J)) END DO FLUX = (FLUXA - FLUXB) / 2 FLUXERROR = RMS / ( 2 * 1.414 ) WRITE (lu,'(9x,''A beam flux = '',1pg12.6,''K'', & 3x,''B beam flux = '',1pg12.6,''K''/ & 9x,''average flux= '',1pg12.6,'' +- _'')') & FLUXA, FLUXB, FLUX * ELSE * single beam flux fit J = 0 DO I = MIDBEAM - 5, MIDBEAM + 5 J = J + 1 CALL GFLUX (ARRAY,I,FLUXESTA(J)) END DO FLUX = FLUXESTA(1) DO J = 2, 11 FLUX = MAX (FLUX,FLUXESTA(J)) END DO FLUXERROR = RMS / 2 WRITE (lu,'(9x,''flux = '',1pg12.6,'' +- _'')') FLUX END IF * WRITE (lu,'(1pg12.6,'' K'')') FLUXERROR IF (FLUX .LT. 5 * FLUXERROR) THEN WRITE (lu,'(9x,''flux less than 5 * flux error :''/ & 3x,''A beam flux estimates (K) :''/ & 3x,6(1pg10.3,2x) / 3x,5(1pg10.3,2x))') & (FLUXESTA(I), I = 1,11) IF (DUALBEAM) THEN WRITE (lu,'(3x,''B beam flux estimates (K) :''/ & 3x,6(1pg10.3,2x) / 3x,5(1pg10.3,2x))') & (FLUXESTB(I), I = 1,11) END IF END IF RETURN END * SUBROUTINE LNOFF (N,DATA,RMS) ************************************************************************ * remove slope between start and end points of data in array * inputs : * N integer*2 last valid point in array * DATA real*4 data array * outputs : * DATA real*4 data array with slope removed * RMS real*4 RMS scatter in first ten and last ten points * IMPLICIT REAL*8 (A-H,O-Z) REAL*4 DATA(N), RMS INTEGER*2 ISTART(2), IEND(2) * * use first ten and last ten points of data for slope * linear regression : MRANGE = 10 SX = 0.0 SXX = 0.0 SY = 0.0 SXY = 0.0 DO K = 1,2 IF (K .EQ. 1) THEN ISTART(K) = 1 IEND(K) = ISTART(K) + MRANGE - 1 ELSE ISTART(K) = N - MRANGE + 1 IEND(K) = N END IF * DO M = ISTART(K), IEND(K) X = M SX = SX + X SXX = SXX + X**2 SY = SY + DATA(M) SXY = SXY + X * DATA(M) END DO END DO * * regression coefficients : M2 = MRANGE * 2 SLOPE = ( SXY - SX*SY/M2 ) / ( SXX - SX*SX/M2 ) CONSTANT = SY/M2 - SLOPE*SX/M2 * subtract linear drift : DO I = 1, N DATA(I) = DATA(I) - ( SLOPE*I + CONSTANT ) END DO * * find rms noise over fitted region : VAR = 0.0 DO K = 1,2 DO M = ISTART(K), IEND(K) VAR = VAR + DATA(M) * DATA(M) END DO END DO RMS = SQRT (VAR/M2) RETURN END * SUBROUTINE GFLUX (DATA,MIDBEAM,FLUX) ****************************************** * calculate flux by fitting a gaussian * inputs : * DATA real*4 data array * MIDBEAM integer*2 index of mid beam = center of Gaussian * output : * FLUX real*4 flux calculated from Gaussian fit * * number of samples per half power beamwidth PARAMETER (NSPBW = 12) * REAL DATA(1) * IRANGE = NSPBW / 2 SFLUX = 0.0 SBMSQ = 0.0 DO I = MIDBEAM-IRANGE, MIDBEAM+IRANGE * offset from center of beam in sampling units OFFSET = FLOAT (I - MIDBEAM) / FLOAT (IRANGE) * normalised beam height at offset - use eqn for Gaussian BEAM = EXP (-0.693 * OFFSET**2) SFLUX = SFLUX + BEAM * DATA(I) SBMSQ = SBMSQ + BEAM**2 END DO * flux in data input units FLUX = SFLUX / SBMSQ RETURN END ************************ INCLUDE DSSUB.FTN INCLUDE HALIM.FTN INCLUDE RDSYS.FTN INCLUDE UPCASE.FTN INCLUDE COMMONLUS.SUB INCLUDE LOSET.SUB INCLUDE OBSOK.SUB ************************