C @(#)refrac.for 17.1.1.1 (ESO-IPG) 01/25/02 17:56:14 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C PROGRAM REFRAC C**************************************************************************** C* * C* Usage: @@ P2:REFRAC infile outfile * C* * C* This program can only be used in MIDAS version may'91 * C* in interactive mode (graphic window must have been opened) * C* * C* * C* Description: * C* * C* This program calculates the atmosperic differential refraction * C* and the light losses of the observed spectra depending on the * C* airmass. The input file is the observed spectra, the output file * C* the calculated spectra which would have been obtained without any * C* light losses (atm.refraction, seeing, slit width). * C* The parallactic angel (slit/elongated aperture perpendicular to horizon* C* will be calculated at any given time. Using this angle at observation * C* ligth losses will be minimized. * C* * C* See also the documentation below the subroutine CALCREC in this file * C* as well as the documentation file REFRAC.DOC * C* * C* VERSION DATE AUTHOR * C* ------------------------------------------------------------------------ C* 1.0 Jan'94 P.W.A. GOERDT * C* * C**************************************************************************** C IMPLICIT NONE C INTEGER ISTAT,NAXISA,NPIXA(2),MODE INTEGER MADRID(1),UNIT,KNUL,ACTS INTEGER*8 PNTRA,PNTRB INTEGER IMNOA,IMNOB,CUNITA C DOUBLE PRECISION STARTA(2),STEPA(2) C REAL WAVE(200000),FAKTOR(200000) C CHARACTER*60 FRAMEA,FRAMEB, CMODE CHARACTER*72 IDENTA CHARACTER*80 STRING C C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/ MADRID COMMON /MOD/ MODE INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C INITIALLISIEREN CALL STSPRO('REFRAC') C C GET NAME OF INPUT FILE CALL STKRDC('IN_A',1,1,60,ACTS,FRAMEA,UNIT,KNUL,ISTAT) CALL STKRDC('IN_B',1,1,60,ACTS,CMODE,UNIT,KNUL,ISTAT) CALL STKRDC('OUT_A',1,1,60,ACTS,FRAMEB,UNIT,KNUL,ISTAT) C C MODE = 0 C TYPE*,'CMODE=',CMODE WRITE(STRING,*) 'CMODE=',CMODE IF (CMODE(1:1) .EQ. 'A') MODE = 1 C C OPEN INPUT IMAGE NAME CALL STIGET (FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,2,NAXISA, & NPIXA,STARTA,STEPA,IDENTA,CUNITA,PNTRA,IMNOA,ISTAT) C C CHECK IMAGE TYPE IF(NAXISA.NE.1) CALL STETER(1,'WRONG IMAGE TYPE!! ') C C CALCULATION OF ATMOSPHERIC REFRACTION CALL CALCREC(STARTA(1),STEPA(1),NPIXA(1),WAVE,FAKTOR) C C CREATE NEW FRAMES CALL STIPUT (FRAMEB,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,1,NPIXA, & STARTA,STEPA,IDENTA,CUNITA,PNTRB,IMNOB,ISTAT) C C FILL FRAMES WITH VALUES MULTIPLIED BY FAKTOR CALL FRAMFILL (MADRID(PNTRA),MADRID(PNTRB),NPIXA,FAKTOR) C C Disconnect from MIDAS enviroment and release keywords CALL STSEPI END SUBROUTINE FRAMFILL(A,B,NPIX,FAKTOR) C REAL A(1),B(1),FAKTOR(*) INTEGER NPIX(2),I C C POINTER CALCULATION DO 200 I=1,NPIX(1) IF (FAKTOR(I).LT.1E-4) THEN CALL STETER(9,'Error: Division by a too small value') ENDIF B(I)=A(I) / FAKTOR(I) 200 CONTINUE C C RETURN END SUBROUTINE CALCREC(START,STEP,NPIX,WAVE,RFLUX) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) EXTERNAL GAUSS ************************************************************************ * * * Written by P.W.Alexander Goerdt, Goettingen (FRG), 1994 * * * * This program calculates the differentiel atmospheric refraction * * depending on the zenith distance of the targets. * * The difference DELTA (in arcsec) of a given wavelength will be * * computed as well as the parallactic angle of the slit (degree), * * which can be rotated to be perpendicular to the horizon * * (North--South = 0.), to get the whole wavelength region of the * * light source into the slit. * * Input data are the slit size (in arcsec) and the astronomical * * angle you have observed the object with and the seeing. * * Hence the flux we will get from a point source observed with a given* * angle will be integrated. * * * * References : * * 1) C.W. Allen, Astrophy. Quantities 3rd edition, p.118ff * * 2) De Ball, Lehrbuch d. sph. Astronomie, Leipzig 1912, p.46 * * 3) Duffett-Smith, Practical astr. with your calculater * * 4) Filippenko, A., P.A.S.P.(1982), Vol. 94, p.715 * ************************************************************************ DOUBLE PRECISION BRECH(2),BLAM,BLHT,BALL,BALLN,DATEJUL,LSTD DOUBLE PRECISION PI,LTR,DECR,RECR,PHR,XLONG,ZENIT DOUBLE PRECISION START,STEP,XLAMBDA REAL ABWEICH(200000) REAL RFLUX(*),WAVE(*) CHARACTER*80 STRING COMMON /MOD/ MODE PI = ACOS(-1.D0) IZAHL = 0 90 CALL STTPUT(' ',ISTAT) CALL STTPUT('Select the unit of your input wavelength!',ISTAT) CALL STTPUT( ' UNIT PRESS ',ISTAT) CALL STTPUT( ' [ANGSTROM] (10**-10m) 1 ',ISTAT) CALL STTPUT( ' [nm] (10** -9m) 2 ',ISTAT) CALL STTPUT( ' [METER] (1m) 3 ',ISTAT) CALL STTPUT( ' Other Unit (?) 4 ',ISTAT) STRING(1:) = 'This program is only valid in the optical '// + 'wavelength region!' CALL GETI(STRING,IUNIT,MODE,1) IF(IUNIT.EQ.1) THEN WUNIT = 1. ELSEIF(IUNIT.EQ.2) THEN WUNIT = 10. ELSEIF(IUNIT.EQ.3) THEN WUNIT = 10.E10 ELSEIF(IUNIT.EQ.4) THEN CALL STTPUT('INPUT ONLY THE EXPONENT OF THE REQUIRED UNIT',ISTAT) CALL GETI(' e.g.: km: PRESS 3 (1km = 10**3m)',IEXPO,MODE,2) WUNIT = 10**(IEXPO+10) ELSE CALL STTPUT(' WRONG INPUT ! TRY AGAIN ',ISTAT) GOTO 90 ENDIF CALL INPUTORT(ALTI,TEMP,WDD,PHG,XLONG) CALL INPUTOBS(RECH,DECG,SIGMA,WAVENORM,XSLIT,YSLIT,ANGLE) CALL INPUTDATE(JAHR,MONAT,NTAG,GMTH) CALL HOEHE(ALTI,P) CALL STDTORAD(RECH,RECR) CALL GRADTORAD(PHG,PHR) CALL DECTORAD(DECG,DECR) CALL STDTODEZ(GMTH,GMTD) CALL GETR('EXPOSURE TIME (min)?',BELT,MODE,1) GMTD = GMTD + BELT / 120. IF(GMTD.GE.24.) THEN GMTD = GMTD - 24. NTAG = NTAG + 1 ENDIF TAG = NTAG + (GMTD / 24.) XLOND = XLONG / 15. CALL JULIENDAY(JAHR,MONAT,TAG,DATEJUL,TAGDJAHRES) CALL SIDEREAL(GMTD,JAHR,XLOND,DATEJUL,TAGDJAHRES,LSTD) CALL DEZTOSTD(LSTD,STH) LTR = LSTD * 2. * PI / 24. CALL PARALLAXE(RECR,DECR,PHR,LTR,ZENIT,WINKEL) CALL AIRMASS(ZENIT,AM) C** CALCULATION OF THE REFRACTION INDEX WNORM = WAVENORM * WUNIT / 1.E4 CALL REFLAM(WNORM,BLAM) CALL REFLHT(BLAM,TEMP,P,BLHT) CALL REFALL(BLHT,WDD,TEMP,WNORM,BALL) BALLN = BALL C** CALCULATION OF THE REFRACTION INDEX AND THE DIFFERENTIEL C** DIFFERENCE IN ARCSEC FOR ALL WAVELENGTH DO 100 NRLAMBDA = 0,(NPIX-1),1 XLAMBDA = START + ( STEP * NRLAMBDA) IZAHL=IZAHL+1 WAVEL = XLAMBDA * WUNIT / 1.E4 IF(WAVEL.LE.0.) THEN CALL STTPUT('ALL WAVELENGTHS SHOULD BE POSITIVE',ISTAT) WRITE(STRING,*) (NRLAMBDA+1),'. WAVELENGTH INPUT: ',XLAMBDA CALL STETER(9,STRING) ENDIF CALL REFLAM(WAVEL,BLAM) CALL REFLHT(BLAM,TEMP,P,BLHT) CALL REFALL(BLHT,WDD,TEMP,WAVEL,BALL) BRECH(1)=BALLN BRECH(2)=BALL CALL DIFFERENZ(BRECH,ZENIT,DELTA) CALL INTEGRAL(XSLIT,YSLIT,DELTA,SIGMA,ANGLE,WINKEL,FLUX) WAVE(IZAHL) = XLAMBDA ABWEICH(IZAHL) = DELTA RFLUX(IZAHL) = FLUX NZAHL = IZAHL 100 CONTINUE CALL OUTPUT(WAVE,ABWEICH,RFLUX,NZAHL,AM,WINKEL,STH,ANGLE) END SUBROUTINE JULIENDAY(JAHR,MONAT,TAG,DATEJUL,TAGDJAHRES) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) DOUBLE PRECISION DATEJUL COMMON /MOD/ MODE JJ = JAHR MM = MONAT IF(MONAT.LE.2) THEN JJ = JAHR - 1 MM = MONAT + 12 ENDIF GREGOR = JAHR + ((MONAT-1)+AINT(TAG)/31)/12 IF(GREGOR.GT.1582.791) THEN A = AINT(1.*JJ/100) B = 2 - A + AINT(A/4) ELSE B = 0 ENDIF C = AINT(365.25*JJ) D = AINT(30.6001*(MM+1)) DAYS=D+AINT(TAG) DATEJUL = B + C + D + TAG + 0.5 DATEJUL = DATEJUL + 1720994. C CALL STTPUT( 'JULIEN DATE = ' , DATEJUL,ISTAT) NSCHALT=0 IF(MOD(JAHR,4).EQ.0.AND.MOD(JAHR,100).NE.0. *OR.MOD(JAHR,400).EQ.0) THEN NSCHALT = 1 ENDIF C CALL STTPUT('IST ',JAHR,' EIN SCHALTJAHR ? (1=JA)',NSCHALT,ISTAT) IF(MONAT.GE.3.AND.MONAT.LE.12) THEN XM = AINT((MONAT + 1) * 30.6) IF(NSCHALT.EQ.0) XM = XM - 63. IF(NSCHALT.EQ.1) XM = XM - 62. ELSEIF(MONAT.LE.2.AND.MONAT.GE.1) THEN XM = (MONAT - 1) IF(NSCHALT.EQ.0) XM = XM * 63. IF(NSCHALT.EQ.1) XM = XM * 62. XM = AINT(XM/2) ELSE CALL STETER(9,'Wrong input for the month!') ENDIF TAGDJAHRES = XM + AINT(TAG) C CALL STTPUT( 'DAY OF THE YEAR : ',TAGDJAHRES,ISTAT) END SUBROUTINE SIDEREAL(GMTD,JAHR,XLOND,DATEJUL,TAGDJAHRES,LSTD) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) DOUBLE PRECISION DATEJB,DATEJUL,A,B,C,E,F,G,T,R,U,GSTD,LSTD PARAMETER(A=0.0657098D0,C=1.002738D0,E=6.6460656D0) PARAMETER(F=2400.051262D0,G=0.00002581D0) COMMON /MOD/ MODE CALL JULIENDAY(JAHR,1,0.,DATEJB,TAG) T = (DATEJB - 2415020.0) / 36525.0 R = E + (F*T) + (G*T**2.D0) U = R - (24.D0*(JAHR-1900.D0)) B = 24.D0 - U GSTD = TAGDJAHRES * A - B + (GMTD * C) IF(GSTD.LT.0.) GSTD = GSTD + 24. IF(GSTD.GT.24.) GSTD = GSTD - 24. C CALL STTPUT( 'GREENWICH SIDERIAL TIME :',GSTD,ISTAT) LSTD = GSTD - XLOND IF(LSTD.LT.0.) LSTD = LSTD + 24. IF(LSTD.GT.24.) LSTD = LSTD - 24. C CALL STTPUT( 'LOCAL SIDERIAL TIME :',LSTD,ISTAT) END SUBROUTINE AIRMASS(ZENIT,AM) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) CHARACTER*80 STRING DOUBLE PRECISION A,B,C,SECZ,ZENIT PARAMETER(A=0.0018167,B=0.002875,C=0.0008083) COMMON /MOD/ MODE SECZ = 1 / ZENIT C CALL STTPUT( 'SEC Z = ',SECZ,ISTAT) AM = SECZ - A*(SECZ-1.) - B*(SECZ-1.)**2 - C*(SECZ-1.)**3 C CALL STTPUT('AIRMASS = ',AM,ISTAT) IF(AM.LT.0.) THEN CALL STTPUT('TARGET NOT OBSERVABLE AT THIS TIME !',ISTAT) WRITE(STRING,456) AM CALL STTPUT(STRING,ISTAT) 456 FORMAT('AIRMASS =',F12.6) CALL STETER(9,'Sorry...') ENDIF END SUBROUTINE DIFFERENZ(BRECH,ZENIT,DELTA) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) DOUBLE PRECISION BRECH(2),ZENIT,PI COMMON /MOD/ MODE PI=ACOS(-1.D0) DO 100 I=1,2 BRECH(I) = BRECH(I) * 1D-6 + 1. 100 CONTINUE ZENRAD = ACOS (ZENIT) DELTA = (BRECH(1) - BRECH(2)) / (BRECH(1) * BRECH(2))*TAN(ZENRAD) ARCSEC = 2. * PI / 360. / 3600. DELTA = DELTA / ARCSEC END SUBROUTINE INPUTORT(ALTI,TEMP,WDD,PHG,XLONG) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) DOUBLE PRECISION XLONG COMMON /MOD/ MODE CHARACTER*9 H 10 FORMAT(A) 111 CALL STTPUT('Which Telescope? ',ISTAT) CALL STTPUT('CALAR ALTO (Spain) : CA ',ISTAT) CALL STTPUT('KITT PEAK (Arizona) : KP ',ISTAT) CALL STTPUT('LA PALMA (Kanary I) : LP ',ISTAT) CALL STTPUT('LA SILLA (Chile) : LS ',ISTAT) CALL STTPUT('Mc DONALD (Texas) : MD ',ISTAT) CALL STTPUT('OTHER OBS. SITE : OS ',ISTAT) CALL GETC('?',H, MODE, 1) IF(INDEX(H,'CA').GT.0.OR.INDEX(H,'ca').GT.0) THEN PHG = 37.23 XLONG = 2.5367 ALTI = 2168. CALL STTPUT(' ',ISTAT) CALL STTPUT('TELESCOPE : CALAR ALTO',ISTAT) ELSEIF(INDEX(H,'LS').GT.0.OR.INDEX(H,'ls').GT.0) THEN PHG =-29.2567 XLONG = 70.73 ALTI = 2347. CALL STTPUT(' ',ISTAT) CALL STTPUT('TELESCOPE : LA SILLA',ISTAT) ELSEIF(INDEX(H,'MD').GT.0.OR.INDEX(H,'md').GT.0) THEN PHG = 30.67167 XLONG = 104.021667 ALTI = 2075. CALL STTPUT(' ',ISTAT) CALL STTPUT('TELESCOPE : Mc DONALD',ISTAT) ELSEIF(INDEX(H,'KP').GT.0.OR.INDEX(H,'kp').GT.0) THEN PHG = 31.96 XLONG = 111.6 ALTI = 2120. CALL STTPUT(' ',ISTAT) CALL STTPUT('TELESCOPE : KITT PEAK',ISTAT) ELSEIF(INDEX(H,'LP').GT.0.OR.INDEX(H,'lp').GT.0) THEN PHG = 28.7594 XLONG = 17.87611 ALTI = 2369. CALL STTPUT(' ',ISTAT) CALL STTPUT('TELESCOPE : LA PALMA',ISTAT) ELSEIF(INDEX(H,'OS').GT.0.OR.INDEX(H,'os').GT.0) THEN CALL STTPUT ('coordinates of observing site: ',ISTAT) CALL GETR('LATITUDE (-90. to 90. DEGREE)', PHG, MODE, 2) RR = SNGL(XLONG) CALL GETR( + 'LONGITUDE (-180 to 180 degree, east of Greenwich negative!)', + RR, MODE, 3) CALL GETR('altitude of observing site [meter] ? ',ALTI, + MODE, 4) ELSE CALL STTPUT(H,ISTAT) CALL STTPUT('Error: Wrong input Telescope, try again:',ISTAT) GOTO 111 ENDIF CALL GETR('AIR TEMPERATURE [C] ? ',TEMP, MODE, 5) CALL DRUCKGRAF CALL GETR('WATER VAPOUR PRESSURE [Torr (mmHg)] ? ',WDD, MODE, 6) END SUBROUTINE INPUTOBS(RECH,DECG,SIGMA,WAVENORM,XSLIT,YSLIT,ANGLE) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) COMMON /MOD/ MODE CALL STTPUT('Coordinates of Target',ISTAT) CALL GETR('RA. (HH.MMSS) ',RECH, MODE, 7) RECH = RECH * 10000. CALL GETR('DECLINATION (GG.MMSS) ',DECG, MODE, 8) DECG = DECG * 10000. CALL GETR('wavelength [in units you selected before] centered' + //' at slit?',WAVENORM, MODE, 9) C WAVENORM=5550 CALL GETR('SEEING [arcsec]?',SEEING, MODE, 10) SIGMA = SEEING / 2.354 CALL GETR('SLIT ANGLE (0-180; 0.& 180.= NORTH--SOUTH) ?',ANGLE, + MODE, 11) CALL GETR('SLIT LENGTH [arcsec]? ',YSLIT, MODE, 12) CALL GETR('SLIT WIDTH [arcsec]?',XSLIT, MODE, 13) END SUBROUTINE INPUTDATE(JAHR,MONAT,NTAG,GMTH) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) COMMON /MOD/ MODE CALL STTPUT('TIME OF OBSERVATION :',ISTAT) CALL GETI('YEAR ? ',JAHR, MODE, 3) CALL GETI(' MONTH ?',MONAT, MODE, 4) CALL GETI('DAY OF MONTH ?',NTAG, MODE, 5) CALL GETR('UNIVERSAL TIME (GMT) [HH.MMSS] ? ',GMTH, MODE, 14) GMTH = GMTH * 10000. END SUBROUTINE REFLAM(WAVEL,BLAM) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) ************************************************************************ * Calculation of the refraction index depending on the wavelength * * for the standard atmosphere * * air pressure : 1013.25 hPa (760 Torr) * * temperature : 15.00 C * * water vapour pressure : 0.00 hPa * ************************************************************************ DOUBLE PRECISION BLAM,TERM1,TERM2,A,B,C,D,E PARAMETER(A=64.328,B=29498.1,C=146,D=255.4,E=41) COMMON /MOD/ MODE TERM1 = A + B / (C - (1. / WAVEL)**2.) TERM2 = D / (E - (1. / WAVEL)**2.) BLAM = TERM1 + TERM2 END SUBROUTINE REFLHT(BLAM,TEMP,P,BLHT) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) ************************************************************************ * Calculation of the refraction index taken into account the * * air pressure and the temperature * ************************************************************************ DOUBLE PRECISION BLAM,BLHT,A,B,C,D,TERM1,TERM2 PARAMETER(A=1.049,B=0.0157,C=720.883,D=0.003661) COMMON /MOD/ MODE TERM1 = P * (1. + (A - B * TEMP) * 1.E-6 * P) TERM2 = C * (1. + D * TEMP) BLHT = BLAM * TERM1 / TERM2 END SUBROUTINE REFALL(BLHT,WDD,TEMP,WAVEL,BALL) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) ************************************************************************ * Calculation of the refraction index taken into account the * * water vapour pressure * ************************************************************************ DOUBLE PRECISION BLHT,BALL,TERM1,TERM2,A,B,C PARAMETER(A=0.0624,B=0.00068,C=0.003661) COMMON /MOD/ MODE TERM1 = A - B / WAVEL**2. TERM2 = 1. + C * TEMP BALL = BLHT - TERM1 / TERM2 * WDD END SUBROUTINE OUTPUT(WAVE,ABWEICH,RFLUX,NZAHL,AM,WINKEL,STH,ANGLE) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) INTEGER NZAHL CHARACTER*80 STRING REAL ABWEICH(*),RFLUX(*),WAVE(*) COMMON /MOD/ MODE 10 FORMAT(A18,F8.3,3X,A18,F8.3) 11 FORMAT(A18,F8.1) 12 FORMAT(A26) 14 FORMAT(A18,F8.4) 19 FORMAT(1X,I5,1X,F9.1,4X,2(F9.5,5X)) DO 100 I=1,2 C CALL STTPUT( 'N(',LAMBDA(I),')=',BRECH(I),ISTAT) 100 CONTINUE STERNZ = STH / 10000. C** TERMINAL OUTPUT WRITE(STRING,14) 'LS TIME (HH.MMSS)=',STERNZ CALL STTPUT(STRING,ISTAT) WRITE(STRING,10) ' AIRMASS =',AM CALL STTPUT(STRING,ISTAT) WRITE(STRING,12) ' -------------------------' CALL STTPUT(STRING,ISTAT) WRITE(STRING,*) ' ' CALL STTPUT(STRING,ISTAT) WRITE(STRING,11) ' ANGLE (DEGREE) = ',WINKEL CALL STTPUT(STRING,ISTAT) WRITE(STRING,12) ' --------------------------' WRITE(STRING,*) ' ' CALL STTPUT(STRING,ISTAT) C** OUTPUT ON FILE UNIT 12 C WRITE(12,10)' BEST ANGLE :',WINKEL,'CHOOSEN ANGLE :', C $ ANGLE C WRITE(12,10)'AIRMASS :',AM C WRITE(12,*) ' NR. LAMBDA DELTA [ARCSEC] FLUSS [%] ' C WRITE(12,*) NZAHL C DO 200 I=1,NZAHL C WRITE(12,19) I,WAVE(I),ABWEICH(I),RFLUX(I) 200 CONTINUE END SUBROUTINE HOEHE(ALTI,P) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) ************************************************************************ * Input : Altitude [meters] * * Output : pressure in Torr (mmHg) * * This subroutine calculates the air pressure P [hectopascal] as * * a function of the altitude h [m] due to the barometric altitude * * formula: * * P(h) = Po * exp (-(Do * g / Po) * h) * * References: Bergmann/Schaefer Vol.1, page 232 * * Here there is: * * Po : pressure at sea level ==> 101325 pascal * * Do : density at sea level ==> 1.2930 kg/m3 * * g : earth acceleration ==> 9.8100 m/s2 * * h : altitude above sea level in meters * ************************************************************************ PARAMETER(DRUCK=101325.,DICHTE=1.293,G=9.81) COMMON /MOD/ MODE C** CALCULATIONS CONST=DICHTE*G/DRUCK PRESSURE = DRUCK / 100. * EXP ( -CONST * ALTI ) P = PRESSURE * 0.75 C CALL STTPUT( 'AIR PRESSURE [Torr] :',P,ISTAT) END SUBROUTINE GRADTORAD(GRAD,RAD) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) CHARACTER*80 STRING DOUBLE PRECISION RAD,PI COMMON /MOD/ MODE PI=ACOS(-1.D0) IF(GRAD.GT.90.OR.GRAD.LT.-90.) THEN WRITE(STRING,*) 'WRONG INPUT OF LATITUDE OF OBSERVATORY! ' CALL STTPUT(STRING,ISTAT) WRITE(STRING,*) 'DEGREE INPUT= ',GRAD CALL STETER(9,STRING) ENDIF RAD = GRAD / 360. * 2. * PI C CALL STTPUT( 'POLHOEHE = ',RAD END SUBROUTINE STDTORAD(HOUR,RAD) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) CHARACTER*80 STRING DOUBLE PRECISION RAD,PI COMMON /MOD/ MODE PI=ACOS(-1.D0) A1 = AINT(HOUR/10000) B1 = A1 / 24. * 360. IF(A1.GE.24.) THEN WRITE(STRING,*) 'WRONG INPUT OF COORDINATES (HOUR)! ' CALL STTPUT(STRING,ISTAT) WRITE(STRING,*) 'HOUR INPUT = ',A1 CALL STETER(9,STRING) ENDIF A2 = MOD(HOUR,10000.) A3 = AINT(A2/100) IF(A3.GE.60.) THEN WRITE(STRING,*) 'WRONG INPUT OF COORDINATES (MIN)! ' CALL STTPUT(STRING,ISTAT) WRITE(STRING,*) 'MIN. INPUT = ',A3 CALL STETER(9,STRING) ENDIF B2 = A3 / 60. * 15. A4 = MOD(A2,100.) IF(A4.GE.60.) THEN WRITE(STRING,*) 'WRONG INPUT OF COORDINATES (SEC)! ' CALL STTPUT(STRING,ISTAT) WRITE(STRING,*) 'SEC. INPUT = ',A4 CALL STETER(9,STRING) ENDIF B3 = A4 / 3600. * 15. DEZ = B1 + B2 + B3 RAD = DEZ / 360. * 2. * PI END SUBROUTINE DECTORAD(GRAH,RAD) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) CHARACTER*80 STRING DOUBLE PRECISION RAD,PI COMMON /MOD/ MODE PI=ACOS(-1.D0) A1 = AINT(GRAH/10000.) IF(A1.GT.90..OR.A1.LT.-90.) THEN WRITE(STRING,*) 'WRONG INPUT OF DECLINATION ! ' CALL STTPUT(STRING,ISTAT) WRITE(STRING,*) 'DEGREE INPUT =',A1 CALL STETER(9,STRING) ENDIF A2 = MOD(GRAH,10000.) A3 = AINT(A2/100) IF(A3.GE.60.) THEN WRITE(STRING,*) 'WRONG INPUT OF DECLINATION (MM) ! ' CALL STTPUT(STRING,ISTAT) WRITE(STRING,*) 'MM INPUT =',A3 CALL STETER(9,STRING) ENDIF B2 = A3 / 60. A4 = MOD(GRAH,100.) IF(A4.GE.60.) THEN WRITE(STRING,*) 'WRONG INPUT OF DECLINATION (SS) ! ' CALL STTPUT(STRING,ISTAT) WRITE(STRING,*) 'SEC INPUT =',A4 CALL STETER(9,STRING) ENDIF B3 = A4 / 3600. DEZ = A1 + B2 + B3 RAD = DEZ / 360. * 2. * PI END SUBROUTINE STDTODEZ(STD,DEZ) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) CHARACTER*80 STRING COMMON /MOD/ MODE A1 = AINT(STD/10000) A2 = MOD(STD,10000.) A2 = AINT (A2/100.) A3 = MOD(STD,100.) IF(A1.GE.24..OR.A2.GE.60..OR.A1.GE.60.) THEN WRITE(STRING,*) 'WRONG INPUT OF COORDINATES (TIMES) ! ' CALL STETER(9,STRING) ENDIF A3 = A3 / 60. A2 = (A2 + A3) / 60. DEZ = A1 + A2 END SUBROUTINE DEZTOSTD(DEZ,STD) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) DOUBLE PRECISION DEZ COMMON /MOD/ MODE A1 = DINT(DEZ) A2 = (DEZ-A1)*60. A3 = MOD((A2),1.)*6./1000. A2 = AINT (A2)/100. STD = (A1 + A2 + A3) * 10000. IF ((A3*10000.).GE.59.5) THEN A3 = 0. IF(A2.EQ.0.59) THEN A2 = 0. IF(A1.EQ.23) THEN A1 = 0. ELSE A1 = A1 + 1. ENDIF ELSE A2 = A2 + 0.01 ENDIF ENDIF STD = (A1 + A2 + A3) * 10000. END SUBROUTINE PARALLAXE(RECR,DECR,PHR,LTR,ZENIT,WINKEL) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) DOUBLE PRECISION LTR,DECR,RECR,PHR,ZENIT,XNENNER,COSWINKEL,PI COMMON /MOD/ MODE PI=ACOS(-1.D0) ZENIT = (SIN(PHR)*SIN(DECR)+COS(PHR)*COS(DECR)*COS(LTR-RECR)) ZAEHLER =(SIN(PHR)*COS(DECR)-COS(PHR)*SIN(DECR)*COS(LTR-RECR)) XNENNER = SQRT(1 - (ZENIT**2.)) COSWINKEL = ZAEHLER / XNENNER COSW = COSWINKEL WINKEL = ACOS (COSW) * 360. / (2. * PI) ABW = REAL(DMIN1(DABS(RECR-LTR),(2.*PI) - DABS(RECR-LTR))) IF (RECR.GT.LTR) THEN IF((RECR+ABW).GT.(2.*PI).AND.(LTR-ABW).LT.0.) THEN C WRITE(STRING,*) 'R > L, but zero in between! ' ELSE C WRITE(STRING,*) 'OTHER ANGLE, Case 1: (R > L) ' WINKEL = 360. - WINKEL ENDIF ENDIF IF(RECR.LT.LTR) THEN IF((LTR+ABW).GT.(2.*PI).AND.(RECR-ABW).LT.0.) THEN C WRITE(STRING,*) 'OTHER ANGLE, Case 2: (R < L)' WINKEL = 360. - WINKEL ENDIF ENDIF C WRITE(STRING,*) 'WINKEL (GRAD) = ',WINKEL END FUNCTION GAUSS(X,SIGMA) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) PI = ACOS (-1.) C IF (ABS(X/SIGMA).GT.12.) THEN C GAUSS = 0. C ELSE GAUSS= 1/(SQRT(2.*PI)*SIGMA) *EXP(-(X**2/2./SIGMA**2)) C ENDIF END SUBROUTINE TRAPZD(GAUSS,SIGMA,A,B,S,N,IT) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) EXTERNAL GAUSS COMMON /MOD/ MODE IF (N.EQ.1) THEN S=0.5*(B-A)*(GAUSS(A,SIGMA)+GAUSS(B,SIGMA)) IT=1 ELSE TNM=IT DEL=(B-A)/TNM X=A+0.5*DEL SUM=0. DO 11 J=1,IT SUM=SUM+GAUSS(X,SIGMA) X=X+DEL 11 CONTINUE S=0.5*(S+(B-A)*SUM/TNM) IT=2*IT ENDIF RETURN END SUBROUTINE POLINT(XA,YA,N,X,Y,DY) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) PARAMETER (NMAX=10) DIMENSION XA(N),YA(N),C(NMAX),D(NMAX) COMMON /MOD/ MODE NS=1 DIF=ABS(X-XA(1)) DO 11 I=1,N DIFT=ABS(X-XA(I)) IF (DIFT.LT.DIF) THEN NS=I DIF=DIFT ENDIF C(I)=YA(I) D(I)=YA(I) 11 CONTINUE Y=YA(NS) NS=NS-1 DO 13 M=1,N-1 DO 12 I=1,N-M HO=XA(I)-X HP=XA(I+M)-X W=C(I+1)-D(I) DEN=HO-HP IF(DEN.EQ.0.)PAUSE DEN=W/DEN D(I)=HP*DEN C(I)=HO*DEN 12 CONTINUE IF (2*NS.LT.N-M)THEN DY=C(NS+1) ELSE DY=D(NS) NS=NS-1 ENDIF Y=Y+DY 13 CONTINUE RETURN END SUBROUTINE INTEGRAL(XSLIT,YSLIT,DELTA,SIGMA,ANGLE,WINKELOPT,FLUX) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) CHARACTER*80 STRING EXTERNAL GAUSS COMMON /MOD/ MODE PI = ACOS (-1.) THETA = ABS (WINKELOPT - ANGLE) IF (THETA.GT.360.) THEN WRITE(STRING,*) 'PROBLEMS WITH ANGLES !!' CALL STTPUT(STRING,ISTAT) WRITE(STRING,*) 'BEST ANGLE : 0 < A < 360 A =',WINKELOPT CALL STTPUT(STRING,ISTAT) WRITE(STRING,*) 'INPUT ANGLE : 0 < B < 180 B =',ANGLE CALL STETER(9,STRING) ENDIF IF (THETA.GE.180.) THETA = THETA - 180. IF (THETA.GT.90.) THETA = 180. - THETA THETA = THETA * 2. * PI / 360. XDELTA = DELTA * SIN (THETA) YDELTA = DELTA * COS (THETA) C** Limits to prevent an under(over-)flow during integration (QROMB) A = XDELTA - XSLIT / 2. IF ((A/SIGMA).GT.12.) A = (12. * SIGMA) IF ((A/SIGMA).LT.-12.) A = -(12. * SIGMA) B = XDELTA + XSLIT / 2. IF ((B/SIGMA).GT.12.) B = (12. * SIGMA) IF ((B/SIGMA).LT.-12.) B = -(12. * SIGMA) IF(ABS(A-B).LT.10.E-6) THEN XFLUX = 0. ELSE CALL QROMB(GAUSS,SIGMA,A,B,XFLUX) ENDIF A = YDELTA - YSLIT / 2. IF ((A/SIGMA).GT.12.) A = (12. * SIGMA) IF ((A/SIGMA).LT.-12.) A = -(12. * SIGMA) B = YDELTA + YSLIT / 2. IF ((B/SIGMA).GT.12.) B = (12. * SIGMA) IF ((B/SIGMA).LT.-12.) B = -(12. * SIGMA) IF(ABS(A-B).LT.10.E-6) THEN YFLUX = 0. ELSE CALL QROMB(GAUSS,SIGMA,A,B,YFLUX) ENDIF FLUX = XFLUX * YFLUX END SUBROUTINE QROMB(GAUSS,SIGMA,A,B,SS) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) EXTERNAL GAUSS PARAMETER(EPS=1.E-6,JMAX=20,JMAXP=JMAX+1,K=5,KM=4) DIMENSION S(JMAXP),H(JMAXP) COMMON /MOD/ MODE H(1)=1. IT = 0 DO 11 J=1,JMAX CALL TRAPZD(GAUSS,SIGMA,A,B,S(J),J,IT) IF (J.GE.K) THEN L=J-KM CALL POLINT(H(L),S(L),K,0.,SS,DSS) IF (ABS(DSS).LT.EPS*ABS(SS)) RETURN ENDIF S(J+1)=S(J) H(J+1)=0.25*H(J) 11 CONTINUE PAUSE 'Too many steps.' END SUBROUTINE DRUCKGRAF IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) CHARACTER*80 STRING COMMON /MOD/ MODE WRITE(STRING,*) CALL STTPUT(STRING,ISTAT) WRITE(STRING,*) + ' Sketch to determine the water vapour pressure' CALL STTPUT(STRING,ISTAT) WRITE(STRING,*) + ' Saturation curve = 100% relative humidity' CALL STTPUT(STRING,ISTAT) WRITE(STRING,16) CALL STTPUT(STRING,ISTAT) WRITE(STRING,1) CALL STTPUT(STRING,ISTAT) WRITE(STRING,2) CALL STTPUT(STRING,ISTAT) WRITE(STRING,3) CALL STTPUT(STRING,ISTAT) WRITE(STRING,4) CALL STTPUT(STRING,ISTAT) WRITE(STRING,5) CALL STTPUT(STRING,ISTAT) WRITE(STRING,6) CALL STTPUT(STRING,ISTAT) WRITE(STRING,7) CALL STTPUT(STRING,ISTAT) WRITE(STRING,8) CALL STTPUT(STRING,ISTAT) WRITE(STRING,9) CALL STTPUT(STRING,ISTAT) WRITE(STRING,10) CALL STTPUT(STRING,ISTAT) WRITE(STRING,11) CALL STTPUT(STRING,ISTAT) WRITE(STRING,12) CALL STTPUT(STRING,ISTAT) WRITE(STRING,13) CALL STTPUT(STRING,ISTAT) WRITE(STRING,14) CALL STTPUT(STRING,ISTAT) WRITE(STRING,15) CALL STTPUT(STRING,ISTAT) WRITE(STRING,17) CALL STTPUT(STRING,ISTAT) WRITE(STRING,18) CALL STTPUT(STRING,ISTAT) WRITE(STRING,19) CALL STTPUT(STRING,ISTAT) 1 FORMAT(5X,'|',71X,'|') 2 FORMAT('P |',71X,'*') 3 FORMAT('R |',68X,'*',2X,'|') 4 FORMAT('E 15-',65X,'*',5X,'-') 5 FORMAT('S |',62X,'*',8X,'|') 6 FORMAT('S |',59X,'*',11X,'|') 7 FORMAT('U |',55X,'*',15X,'|') 8 FORMAT('R 10-',51X,'*',19X,'-') 9 FORMAT('E |',46X,'*',24X,'|') 10 FORMAT(5X,'|',40X,'*',30X,'|') 11 FORMAT('T |',34X,'*',36X,'|') 12 FORMAT('O 5 -',27X,'*',43X,'-') 13 FORMAT('R |',19X,'*',51X,'|') 14 FORMAT('R |',10X,'*',60X,'|') 15 FORMAT(5X,'*',71X,'|') 16 FORMAT(2X,'20 -',6(11('-'),'|')) 17 FORMAT(3X,'0 |',6(11('-'),'|')) 18 FORMAT(4X,'-10',9X,'-5',11X,'0',11X,'5',10X,'10',10X,'15', + 10X,'20') 19 FORMAT(34X,'TEMPERATURE [C]') WRITE(STRING,*) + ' e.g.: at 70% humidity and T=10 C the pressure would' + //' be about 6 Torr' CALL STTPUT(STRING,ISTAT) END SUBROUTINE GETR(PROMPT, VALUE, MODE, RANK) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) CHARACTER*(*) PROMPT REAL VALUE INTEGER MODE, RANK INTEGER ACTS, KUN, KNUL, STAT IF (MODE .EQ. 0) THEN CALL STKPRR(PROMPT,'INPUTR',1,1,ACTS,VALUE,KUN,KNUL,STAT) ELSE CALL STKRDR('INPUTR',RANK,1,ACTS,VALUE,KUN,KNUL,STAT) ENDIF END SUBROUTINE GETI(PROMPT, VALUE, MODE, RANK) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) CHARACTER*(*) PROMPT INTEGER VALUE ,MODE, RANK INTEGER ACTS, KUN, KNUL, STAT IF (MODE .EQ. 0) THEN CALL STKPRI(PROMPT,'INPUTI',1,1,ACTS,VALUE,KUN,KNUL,STAT) ELSE CALL STKRDI('INPUTI',RANK,1,ACTS,VALUE,KUN,KNUL,STAT) ENDIF END SUBROUTINE GETC(PROMPT, VALUE, MODE, RANK) IMPLICIT INTEGER(I-N) IMPLICIT REAL(A-H,O-Z) CHARACTER*(*) PROMPT CHARACTER*9 VALUE INTEGER MODE, RANK INTEGER ACTS, KUN, KNUL, STAT IF (MODE .EQ. 0) THEN CALL STKPRC(PROMPT,'INPUTC',1,1,9,ACTS,VALUE,KUN,KNUL,STAT) ELSE CALL STKRDC('OUT_B',1,RANK,9,ACTS,VALUE,KUN,KNUL,STAT) ENDIF END