SUBROUTINE QCAL( IBASE, HA, DEC, WAVE, ANGDIAM, VISEST ) C----------------------------------------------------------------------- C Calculates the visibility squared, given and angular diameter, C wavelength and baseline components. The input diameter is the C uniform disk diameter at 800 nm. C Derived from OICAL, 8/88, D. J. Hutter C Major rewrite, 18 Sept 1988 (D. Mozurkewich) C Changed to v squared 16 Jan 1989 (d.m.) C Added wavelength dependent diameters 28 sep 1990 (dm) C----------------------------------------------------------------------- C inputs C IBASE - Baseline number, i.e. 1 = NS, 2 = ES, >3 => VB C HA - Hour Angle in hours. C DEC - Declination in degrees. C WAVE - Wavelength in nm. C ANGDIAM - Estimated angular diameter in mas C outputs C VISEST - Expected visibility squared C C The baseline components are stored in file X:\baseline.dat C To add a baseline, add it to the file and increase MAXB C----------------------------------------------------------------------- IMPLICIT UNDEFINED(A-Z) SAVE INTEGER *4 MAXB REAL *4 PI, RAD, LON PARAMETER (MAXB= 39) PARAMETER (PI=3.141592653589793) PARAMETER (RAD = 180./PI ) PARAMETER (LON=118.059167) C INTEGER *4 IBASE, IERR, I, J REAL *4 HA, DEC, WAVE, ANGDIAM, VISEST, TEMP, SCALE C INTEGER *4 IB REAL*4 BASEX(MAXB), BASEY(MAXB), BASEZ(MAXB) REAL*4 ARG, BJ1, TIHA, U, V, R, COSDEC, SINDEC REAL*4 BSXI, BSYI, STXROT, STYROT, BXY, BZ, BPROJ CHARACTER *5 TITLE LOGICAL FIRST DATA FIRST / .TRUE. / C----------------------------------------------------------------------- C On the first call, Read in baseline components. C The x, y, z, components are in meters, relative to the south C siderostat and aligned with the earth's rotation at Mount Wilson. C----------------------------------------------------------------------- IF ( FIRST ) THEN FIRST = .FALSE. OPEN (UNIT=56, FILE='X:\BASELINE.DAT', STATUS='OLD', IOSTAT=IERR) IF ( IERR .NE. 0 ) THEN WRITE(6,*) ' ******ERROR****** baseline file is missing ' WRITE(6,*) ' cannot find file X:\baseline.dat' WRITE(6,*) ' ****************************************** ' STOP END IF READ (56,'(/A)',IOSTAT=IERR) TITLE IF ( IERR .NE. 0 ) THEN WRITE(6,*) ' ******ERROR****** baseline file is missing ' WRITE(6,*) ' cannot read file X:\baseline.dat ' WRITE(6,*) ' ****************************************** ' STOP END IF DO 20 I = 1, MAXB READ(56,'(I5,3F15.0)',IOSTAT=IERR) $ J, BASEX(I), BASEY(I), BASEZ(I) IF ( IERR .NE. 0 ) THEN WRITE(6,*) ' ******ERROR****** baseline file is bad ' WRITE(6,*) ' file X:\baseline.dat is too short ' WRITE(6,*) ' ************************************** ' STOP END IF IF ( J .NE. I ) THEN WRITE(6,*) ' BASELINE NUMBERS ARE NOT PROPERLY ', $ 'SEQUENCED. PROBABLY CORRUPT' END IF 20 CONTINUE CLOSE(UNIT=56) C----------------------------------------------------------------------- C Rotate right handed local baseline XYZ to the C left handed frame fixed at 0 longitude C----------------------------------------------------------------------- DO 30 IB=1,MAXB TEMP = BASEX(IB)*COS(LON/RAD) + BASEY(IB)*SIN(LON/RAD) BASEY(IB) = BASEX(IB)*SIN(LON/RAD) - BASEY(IB)*COS(LON/RAD) BASEX(IB) = TEMP 30 CONTINUE END IF C----------------------------------------------------------------------- C End of initiallization C----------------------------------------------------------------------- C Baseline components are in meters, and have to be converted C to millions of wavelengths. C Wavelength is in nanometers. C----------------------------------------------------------------------- IF ( IBASE .GT. MAXB ) THEN WRITE(6,*) ' BASELINE NUMBER ', IBASE, ' IS TOO LARGE! ' WRITE(6,*) ' STAR ASSUMED TO BE UNRESOLVED ' VISEST = 1.0 ELSE C BPROJ = SQRT ( BASEX(IBASE)*BASEX(IBASE) + C $ BASEY(IBASE)*BASEY(IBASE) + C $ BASEZ(IBASE)*BASEZ(IBASE) ) COSDEC = COS(DEC/RAD) SINDEC = SIN(DEC/RAD) BSXI = 1000. * BASEX(IBASE) / WAVE BSYI = 1000. * BASEY(IBASE) / WAVE BZ = 1000. * BASEZ(IBASE) / WAVE BXY = SQRT( BSXI*BSXI + BSYI*BSYI ) C----------------------------------------------------------------------- C Determine (u,v) values. C----------------------------------------------------------------------- TIHA = (15.*HA+LON-90.)/RAD - ATAN2(BSYI,BSXI) U = BXY * COS(TIHA) V = BZ * COSDEC + BXY * SINDEC * SIN(TIHA) R = SQRT ( U*U + V*V ) C----------------------------------------------------------------------- C Calculate the visibility based on the uniform disk diameter. C Note that we take the absolute value of the bessel function to C keep the visibility positive beyond the first minimum. C C Scale the uniform disk diameter from 800 nm to the wavelength C of interest. C----------------------------------------------------------------------- IF ( ABS(WAVE-800.) .LT. 10. ) THEN SCALE = 1.000 ELSE IF ( ABS(WAVE-700.) .LT. 10. ) THEN SCALE = 1.020 ELSE IF ( ABS(WAVE-550.) .LT. 10. ) THEN SCALE = 1.040 ELSE IF ( ABS(WAVE-500.) .LT. 10. ) THEN SCALE = 1.055 ELSE IF ( ABS(WAVE-450.) .LT. 10. ) THEN SCALE = 1.070 ELSE SCALE = 0.91 + 72./WAVE C WRITE(6,*) ' BAD WAVELENGTH FOR SCALING ' END IF C WRITE(6,*) ' ANG DIAMETER SCALING = ', SCALE ARG = 0.015231 * R * ANGDIAM / SCALE IF (ARG .GE. 0.00001 ) THEN CALL BESSJ1(ARG,BJ1) VISEST = 2.0 * ABS(BJ1) / ARG ELSE VISEST = 1.0 END IF VISEST = VISEST * VISEST END IF RETURN END SUBROUTINE BESSJ1(X,BJ1) C----------------------------------------------------------------------- C D.J.Hutter (10/87) C RETURNS THE BESSEL FUNCTION J1(X) FOR ANY REAL X. C TAKEN FROM p.174 OF "NUMERICAL RECIPIES" C----------------------------------------------------------------------- C IMPLICIT UNDEFINED (A-Z) SAVE REAL *4 X, BJ1, AX, Z, XX REAL*4 Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6, * S1,S2,S3,S4,S5,S6 DATA R1, R2, R3, R4, R5, R6 / * 7.2362614232E10, -7895059235., 242396853.1, * -2972611.439, 15704.48260, -30.16036606 / DATA S1, S2, S3, S4, S5, S6 / * 14.4725228442E10, 2300535178., 18583304.74, * 99447.43394, 376.9991397, 1. / DATA P1,P2,P3,P4,P5/1.,0.00183105,-.3516396496E-4, * .2457520174E-5,-.240337019E-6/, * Q1,Q2,Q3,Q4,Q5/.04687499995,-.2002690873E-3, * .8449199096E-5,-.88228987E-6,.105787412E-6/ IF(ABS(X).LT.8.) THEN Y=X**2 BJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))))) * /(S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))))) ELSE AX=ABS(X) Z=8./AX Y=Z**2 XX=AX-2.356194491 BJ1=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y* * (P4+Y*P5))))-Z*SIN(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5))))) * *SIGN(1.,X) END IF RETURN END