C C TRANFER CURVE C Dan Bruton (astro@tamu.edu) C March 27, 1996 C C This program converts a true image of the Sun to an apparent C image for a given temperature profile (input file) using a parabolic C ray path, multilayer atmosphere model. C A portable pixmap (ppm) image file is generated. C C ALT = altitude of the object in arcminutes C AZM = relative azimuth of the object in arcminutes C ALFA = apparent altitude in arcminutes C ALFT = true altitude in arcminutes C RI = refractive index C PRS = pressure in Pascals C TMP = absolute temperature C RC = radius of curvature of the ray in meters C HOBS = height of the observer above the surface of the Earth C ROBJ = radius of the object in arcminutes C RE = radius of the Earth in meters C CV = pixel color RGB values C IMAGE = make image (true or false: 1 or 0) C ITRAN = transfer function (computed or standard: 1 or 0) C BETA = incident angle relative to each layer in radians C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CV(600,600,3) DIMENSION ELV(0:2000), TMP(0:2000), RCV(81,0:2000) DIMENSION ALFA(81,500), ALFT(81,500) DIMENSION EDIS(81), IDUCT(500), AIRM(500) DIMENSION ALT(5), AZM(5) DATA ALT/-20.,-30,-40.,-50.,-60./ DATA AZM/19.,51.,83.,115.,147./ CHARACTER *70 TXT C C PPM IMAGE CONSTANTS AND OTHER INPUT PARAMETERS C HEIGHT, WIDTH, PIXEL DEPTH, PIXELS PER ARCMINUTE C RE=6378140.D0 PI=3.141592654D0 CONV=PI/(180.D0*60.D0) ROBJ=15.D0 N=300 M=500 KMIN=1 KMAX=81 KSTEP=40 GAM=1.D0 MAX=255 RMAX=REAL(MAX) HDEG=1.5D0 PPAM=REAL(N)/(HDEG*60.D0) DELTA=1000.D0 XMAX=30000.D0 HOBS=10.01D0 DELALT=-10.D0 RMX=0.D0 TSUN=5800.D0 ISTND=1 ITRAN=1 IMAGE=1 IPLOT=0 C C READ INPUT PARAMETERS C OPEN(UNIT=20,FILE='param.dat',STATUS='UNKNOWN') READ(20,*) HOBS,ISTND,ITRAN,IMAGE,IPLOT,DELALT,XMAX READ(20,*) TXT READ(20,*) N,M,KMIN,KMAX,KSTEP M=INT(N*6/3) PPAM=REAL(N)/(HDEG*60.D0) CLOSE (20) C C *** MAKE PLOT PARAMETERS *** C IF (IPLOT.EQ.1) THEN N=10 PPAM=0.8333D0 ENDIF C C READ TEMPERATURE PROFILE C PRINT *, ' Reading temperature profile...' READ(*,*) NNN NNN=NNN-1 DO II=0,NNN READ(*,*) ELV(II),TMP(II) c PRINT *, NNN, ELV(II),TMP(II) ENDDO PRINT *, ' Computing radius of curvatures...' AMZ=0.D0 DO II=0,NNN DO K=1,81 WL=380.D0+REAL(K-1)*5.D0 CALL RADCUR(II,RC,AM,ELV,TMP,WL) RCV(K,II)=RC EDIS(K)=0.D0 ENDDO AMZ=AMZ+AM ENDDO PRINT *, ' Vertical Airmass (kg/m^2)',AMZ C C FIND ALFT IN TERMS OF ALFA C PRINT *, ' Finding transfer curve...' DO K=KMIN,KMAX,KSTEP KSTOP=K OPEN(UNIT=20,FILE='test.dat',STATUS='UNKNOWN') 1 FORMAT(F15.5,1X,F15.5) IF (ITRAN.EQ.1) THEN WL=380.D0+REAL(K-1)*5.D0 PRINT *, ' Wavelength (nm) ',WL C FIND THE STARTING LAYER DO III=0,NNN IF (HOBS.GT.ELV(III)) IIS=III+1 ENDDO C RAY TRACE THROUGH THE FIRST LAYER DO J=1,N IF (IPLOT.EQ.1) PRINT *, J II=IIS IDUCT(J)=0 AIRM(J)=0.D0 ALFA(K,J)=(REAL(N/2 - J)/PPAM) BETA=ALFA(K,J)*CONV RC=RCV(K,II)*DCOS(BETA)*DCOS(BETA) X=0.D0 C INITIAL RAY - LEHN'S EQUATION 12 & 17 IF (BETA.GE.0.D0) THEN DIS=DTAN(BETA)**2+2.D0*(ELV(II)-HOBS)*(RC-RE)/ * (RC*RE) IF (DIS.LT.0D0) THEN DIS=DTAN(BETA)**2-2.D0*(HOBS-ELV(II-1))*(RC-RE)/ * (RC*RE) UP=(-DTAN(BETA)-DSQRT(DIS))*(RE*RC/(RC-RE)) IFLAG=0 ELSE UP=(-DTAN(BETA)+DSQRT(DIS))*(RE*RC/(RC-RE)) IFLAG=1 ENDIF ELSE DIS=DTAN(BETA)**2-2.D0*(HOBS-ELV(II-1))*(RC-RE)/ * (RC*RE) IF (DIS.LT.0D0) THEN DIS=DTAN(BETA)**2+2.D0*(ELV(II)-HOBS)*(RC-RE)/ * (RC*RE) UP=(-DTAN(BETA)+DSQRT(DIS))*(RE*RC/(RC-RE)) IFLAG=1 ELSE UP=(-DTAN(BETA)-DSQRT(DIS))*(RE*RC/(RC-RE)) IFLAG=0 ENDIF ENDIF C CONTINUE RAY TRACING THROUGH THE LAYERS THETA=0.D0 H=HOBS DO WHILE ((II.GE.1).AND.(II.LE.NNN).AND.(X.LT.XMAX)) C AIRMASS CALCULATION IF (IFLAG.EQ.1) THEN DZ=ELV(II)-UP**2/(2.D0*RE)-H ELSE DZ=ELV(II-1)-UP**2/(2.D0*RE)-H ENDIF DD=DSQRT(UP*UP + DZ*DZ) AMGAM=2.D0*DASIN(DD/(2.D0*RC)) ADEN=(DEXP(-ELV(II)/8400.D0)+DEXP(-ELV(II-1)/8400.D0))/2.D0 ADEN=1.2250D0*ADEN AIRM(J)=AIRM(J)+ADEN*RC*AMGAM IF (IPLOT.EQ.1) WRITE (20,1) X,H XD=0.D0 HD=H C FILL IN THE GAPS DURING RAY TRACING DO WHILE ((IPLOT.EQ.1).AND.(XD.LT.UP).AND. * ((X+XD).LT.XMAX)) IF (IPLOT.EQ.1) WRITE(20,1) X+XD,HD XD=XD+DELTA HD=((1.D0/RE)-(1.D0/RC))*(XD**2.D0)/2.D0+XD*DTAN(BETA)+H ENDDO C FIND NEW POSITION ANGLE THETA WRT EARTH'S CENTER R1=RE+ELV(II-1) R2=RE+ELV(II) Z=-0.5D0*DCOS(BETA)*UP**2/RC+UP*DTAN(BETA) ARGT=(R1**2 + R2**2 - UP**2 - Z**2)/(2.D0*R1*R2) IF (ARGT.LT.1.D0) THETA=THETA+DACOS(ARGT) C DETERMINE THE NEW PARAMETERS BETA=DATAN(DTAN(BETA)-(UP/RC)+(UP/RE)) ALFT(K,J)=(BETA-THETA)/CONV IF (IFLAG.EQ.1) THEN H=ELV(II) II=II+1 ELSE H=ELV(II-1) II=II-1 ENDIF RC=RCV(K,II)*DCOS(BETA)*DCOS(BETA) X=X+UP C CALCULATE NEW UP VALUE IF ((II.GE.1).AND.(II.LE.NNN)) THEN IF (IFLAG.EQ.1) THEN C LEHN'S EQUATION 12 & 13 DIS=DTAN(BETA)**2+2.D0*(ELV(II)-ELV(II-1))*(RC-RE)/ * (RC*RE) IF (DIS.LT.0D0) THEN UP=2.D0*RE*RC*DTAN(BETA)/(RE-RC) IFLAG=0 ELSE UP=(-DTAN(BETA)+DSQRT(DIS))*(RE*RC/(RC-RE)) IFLAG=1 ZL=ELV(II) ZR=ELV(II-1) ENDIF ELSE C LEHN'S EQUATION 17 & 13 DIS=DTAN(BETA)**2-2.D0*(ELV(II)-ELV(II-1))*(RC-RE)/ * (RC*RE) IF (DIS.LT.0D0) THEN UP=2.D0*RE*RC*DTAN(BETA)/(RE-RC) IFLAG=1 ELSE UP=(-DTAN(BETA)-DSQRT(DIS))*(RE*RC/(RC-RE)) IFLAG=0 ENDIF ENDIF ENDIF ENDDO IF (II.EQ.0) ALFT(K,J)=-1100.D0 IF (X.GE.XMAX) IDUCT(J)=1 ENDDO C C TRANSFER FUNCTION FOR A STANDARD ATMOSPHERE C ELSEIF (ITRAN.EQ.0) THEN DO J=1,N ALFA(K,J)=(REAL(N/2 - J)/PPAM) BETA=(REAL(N/2 - J)/PPAM) BETAD=BETA/60.D0 PR = 99975.D0 TR = 273.D0 TOP=PR*(1.594D-3+(1.96D-4*BETAD)+(2.0D-7*BETAD*BETAD)) BOT=TR*(1.D0+(5.05D-1*BETAD)+(8.45D-2*BETAD*BETAD)) DAB = 60.D0*TOP/BOT ALFT(K,J) = BETA - DAB ENDDO ENDIF C C ADD STADARD ATMOSPHERE DEFLECTION C IF (ISTND.EQ.1) THEN DO J=1,N BETA=ALFT(K,J) BETAD=BETA/60.D0 PR = 99975.D0 TR = 273.D0 TOP=PR*(1.594D-3+(1.96D-4*BETAD)+(2.0D-7*BETAD*BETAD)) BOT=TR*(1.D0+(5.05D-1*BETAD)+(8.45D-2*BETAD*BETAD)) DAB = 60.D0*TOP/BOT ALFT(K,J) = BETA - DAB ENDDO ENDIF ENDDO CLOSE (20) C C READ OR WRITE TRANSFER CURVE TO A FILE C IF (ITRAN.EQ.2) THEN PRINT *, ' Reading transfer curve...' OPEN(UNIT=20,FILE='tc.dat',STATUS='UNKNOWN') READ (20,*) KMIN,KMAX,KSTEP READ (20,*) N, PPAM DO K=KMIN,KMAX,KSTEP DO J=1,N READ(20,1) ALFA(K,J),ALFT(K,J) ENDDO ENDDO CLOSE (20) ELSE PRINT *, ' Writing transfer curve...' OPEN(UNIT=20,FILE='tc.dat',STATUS='UNKNOWN') WRITE (20,*) KMIN,KMAX,KSTEP WRITE (20,*) N, PPAM DO K=KMIN,KMAX,KSTEP DO J=1,N IF (ALFT(K,J).GT.-500.D0) WRITE(20,1) ALFA(K,J),ALFT(K,J) ENDDO ENDDO CLOSE (20) ENDIF C C COLOR BACKGROUND C BSOL=30.D0 c CALL RSBG(CV,BSOL,N,M,PPAM,KMIN,KMAX,KSTEP) C C MAKE IMAGE OF THE OBJECT C USE LIMB DARKENING C IF (IMAGE.EQ.1) THEN PRINT *, ' Making Image...' DO J=1,N IPER=INT(J*10.D0/N) IF(IPER.EQ.(J*10.D0/N)) PRINT *, IPER*10,'%' DO I=1,M DO K=KMIN,KMAX,KSTEP EDIS(K)=0.D0 DO KA=1,5 Y=ALFT(K,J)-ALT(KA)-DELALT X=(REAL(I)/PPAM)-AZM(KA) RHO=DSQRT(X*X+Y*Y) IF ((RHO.LT.ROBJ).AND.(ALFT(K,J).GT.-900.D0)) THEN WL=380.D0+REAL(K-1)*5.D0 CON=1240./8.617D-5 C BLACKBODY RADIATION EDIS1=3.74183D-16*(1./WL**5)/(DEXP(CON/(WL*TSUN))-1.) DARK=1.D0-0.6D0*(1.D0-DSQRT(1.D0-(RHO/ROBJ)**2)) EDIS1=EDIS1*DARK/(4.D0*PI) c *PI*((0.25D0*PI/180.D0)**2) ELSE EDIS1=0.D0 ENDIF TAUN=(283./WL)**4 EDIS(K)=EDIS(K)+EDIS1*DEXP(-TAUN*AIRM(J)/AMZ) ENDDO ENDDO CALL EXYZ(XC,YC,ZC,EDIS) CALL XYZTORGB(XC,YC,ZC,R,G,B) CV(I,J,1)=CV(I,J,1)+R CV(I,J,2)=CV(I,J,2)+G CV(I,J,3)=CV(I,J,3)+B IF (DMAX1(CV(I,J,1),CV(I,J,2),CV(I,J,3)).GT.RMX) THEN RMX=DMAX1(CV(I,J,1),CV(I,J,2),CV(I,J,3)) ENDIF ENDDO ENDDO C C WRITE TO PPM FILE C PRINT *, ' Writing to PPM file...' OPEN(UNIT=20,FILE='temp.ppm',STATUS='UNKNOWN') 1001 FORMAT (A2) WRITE(20,1001) 'P3' 1002 FORMAT (A12) WRITE(20,1002) '# temp.ppm' 1003 FORMAT (I3,1X,I3) WRITE(20,1003) M,N 1004 FORMAT (I3) WRITE(20,1004) MAX DO J=1,N IPER=INT(J*10.D0/N) IF(IPER.EQ.(J*10.D0/N)) PRINT *, IPER*10,'%' DO I=1,M 1005 FORMAT (I3,1X,I3,1X,I3) IR=INT(REAL(MAX)*(255./183.)*(CV(I,J,1)/RMX)**GAM) IG=INT(REAL(MAX)*(CV(I,J,2)/RMX)**GAM) IB=INT(REAL(MAX)*(255./246.)*(CV(I,J,3)/RMX)**GAM) C DRAW ASTRONOMICAL HORIZON IF ((J.EQ.INT(N/2)).AND.(I.LT.10)) THEN IR=255 IG=255 IB=255 ENDIF C DRAW SURFACE OF THE EARTH IF (ALFT(KSTOP,J).LT.-1050.D0) THEN IR=0 IG=0 IB=100 ENDIF C SHOW DUCTING REGION IF (IDUCT(J).EQ.1) IB=100 WRITE(20,1005) IR,IG,IB ENDDO ENDDO CLOSE (20) ENDIF PRINT *, ' Program complete.' END C C********************************************************************* C SUBROUTINE RADCUR(II,RC,AM,ELV,TMP,WL) C C This subroutine calculates the radius of curvature C RC for horizontal rays in the layer II as well as C vertical air mass. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION ELV(0:1000), TMP(0:1000) C IF (II.GT.0) THEN C DRY AIR AND WATER VAPOR PARTIAL PRESSURES PRS=101325.D0 PRW=0.04D0*PRS c TAVE=(TMP(II)+TMP(II-1))/2.D0 c TRG=(TMP(II)-TMP(II-1))/(ELV(II)-ELV(II-1)) c EAVE=(ELV(II)+ELV(II-1))/2.D0 c CON1=9.81D0*28.96D-3/8.314510 c PAVE=PRS*DEXP(-CON1*EAVE/TAVE) c CON2=2.26D-4*28.96D-3/8.314510 c RI=1.D0+CON2*PAVE/TAVE c RKAPPA=(CON2*PAVE/(RI*TAVE*TAVE))*(CON1+TRG) TERM2=0.6839397D0/(130.0D-6 - 1.D0/WL**2) TERM3=0.0045473D0/(38.9D-6 - 1.D0/WL**2) TERM5=6487.31+58.058D6/WL**2-0.71150D12/WL**4+0.08851D18/WL**6 P1=PRS*DEXP(-ELV(II-1)/8400.D0) TERM4=57.90D-10 - 9.3250D-6/TMP(II-1) + 0.25844D-2/TMP(II-1)**2 DS=(P1/TMP(II-1))*(1.D0+P1*TERM4) PW1=PRW*DEXP(-ELV(II-1)/8400.D0) TERM6=1.D0 + 3.7D-6 * PW1 TERM7=-2.37321D-5 + 2.23366D-2/TMP(II-1) TERM8=-710.792D-2/TMP(II-1)**2 + 7.7514E2/TMP(II-1)**3 DW=(PW1/TMP(II-1))*(1.D0+PW1*TERM6*(TERM7+TERM8)) RI1=1.D0+1.0D-10*((2371.34D0+TERM2+TERM3)*DS + TERM5*DW) P2=PRS*DEXP(-ELV(II)/8400.D0) TERM4=57.90D-10 - 9.3250D-6/TMP(II) + 0.25844D-2/TMP(II)**2 DS=(P2/TMP(II))*(1.D0+P2*TERM4) PW2=PRW*DEXP(-ELV(II)/8400.D0) TERM6=1.D0 + 3.7D-6 * PW2 TERM7=-2.37321D-5 + 2.23366D-2/TMP(II) TERM8=-710.792D-2/TMP(II)**2 + 7.7514E2/TMP(II)**3 DW=(PW2/TMP(II))*(1.D0+PW2*TERM6*(TERM7+TERM8)) RI2=1.D0+1.0D-10*((2371.34D0+TERM2+TERM3)*DS + TERM5*DW) DELV=ELV(II)-ELV(II-1) RKAPPA=-2.D0*(RI2-RI1)/(DELV*(RI2+RI1)) RC=1.D0/RKAPPA AM=(DEXP(-ELV(II)/8400.D0)+DEXP(-ELV(II-1)/8400.D0))/2.D0 AM=1.2250D0*AM*DELV C C CIDDOR INDEX OF REFRACTION C XCPPM=450.D0 C INDEX OF REFRACTION FOR DRY AIR SIGMA=1.D0/(WL*1000.D0) RIAS=1.D0+1.D-8*(5792105.D0/(238.0185D0-SIGMA**2) * +167917.D0/(57.362D0-SIGMA**2)) RIAXS=1.D0+(RIAS-1.D0)*(1.D0+0.534D-6*(XCPPM-450.D0)) RIWS=1.D0+1.D-8*1.022D0*(295.235D0+2.6422D0*SIGMA**2 * -0.032380*SIMGA**4 +0.004028*SIGMA**6) C BOTTOM OF THE LAYER CALL DENSITY(DAXS,DW,288.15D0,101325.D0,0.D0,XCPPM,0) CALL DENSITY(DA,DWS,293.15D0,1333.D0,1333.D0,XCPPM,1) P=PRS*DEXP(-ELV(II-1)/8400.D0) T=TMP(II-1) PW=PRW*DEXP(-ELV(II-1)/8400.D0) CALL DENSITY(DA,DW,T,P,PW,XCPPM,2) RI1=1.D0+(DA/DAXS)*(RIAXS-1.D0)+(DW/DWS)*(RIWS-1.D0) C TOP OF THE LAYER P=PRS*DEXP(-ELV(II)/8400.D0) T=TMP(II) PW=PRW*DEXP(-ELV(II)/8400.D0) CALL DENSITY(DA,DW,T,P,PW,XCPPM,2) RI2=1.D0+(DA/DAXS)*(RIAXS-1.D0)+(DW/DWS)*(RIWS-1.D0) C RADIUS OF CURVATURE DELV=ELV(II)-ELV(II-1) RKAPPA=-2.D0*(RI2-RI1)/(DELV*(RI2+RI1)) RC=1.D0/RKAPPA C AIRMASS AM=(DEXP(-ELV(II)/8400.D0)+DEXP(-ELV(II-1)/8400.D0))/2.D0 AM=1.2250D0*AM*DELV ELSE RC=1.0D+20 AM=0.D0 ENDIF RETURN END C C********************************************************************* C C SUBROUTINE DENSITY(DA,DW,T,P,PW,XCPPM,IXW) C C This subroutine calculates the density of moist air C as described by P.E. Ciddor, Applide Optics, Vol. 35, C No. 9., March 20, 1996. MKS units. C C DA = density of dry air component C DW = density of water vapor compoment C T = temperature in Kelvin C P = total pressure C PW = partial pressure of water vapor C XCPPM = CO2 content in parts per million C Fractional Humidity = PW/SVP C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C SATURATION VAPOR PRESSURE OF WATER VAPOR IN AIR C A=1.2378847D-5 B=-1.9121316D-2 C=33.93711047D0 D=-6.3431645D3 SVP=DEXP(A*T*T+B*T+C+D/T) C C ENHANCEMENT FACTOR OF WATER VAPOR IN AIR C ALPHA=1.00062D0 BETA=3.14D-8 GAMMA=5.6D-7 F=ALPHA+BETA*P+GAMMA*(T-273.15D0) C C MOLAR FRACTION OF WATER VAPOR IN MOIST AIR C IF (IXW.EQ.0) THEN XW=0.D0 ELSEIF (IXW.EQ.1) THEN XW=1.D0 ELSE XW=F*PW/P ENDIF C C MOLAR MASS OF DRY AIR AND MOLAR MASS OF WATER VAPOR C RMA=1.D-3*(28.9635D0+12.011D-6*(XC-400.D0)) RMW=0.018015D0 C C COMPRESSIBILITIES OF DRY AIR AND PURE WATER VAPOR C CALL COMPRESS(Z,T,P,XW) C C DENSITY OF MOIST AIR AND ITS COMPONENTS C R=8.314510D0 DA=P*RMA*(1.D0-XW)/(Z*R*T) DW=P*RMW*XW/(Z*R*T) RETURN END C C********************************************************************* C C SUBROUTINE COMPRESS(COMP,T,P,XW) C C This subroutine calculates the compressibility of moist air C as described by P.E. Ciddor, Applide Optics, Vol. 35, C No. 9., March 20, 1996. MKS units. C C T = temperature in Kelvin C P = total pressure C XW = molar fraction of water vapor C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C SATURATION VAPOR PRESSURE OF WATER VAPOR IN AIR C A=1.58123D-6 A1=-2.9331D-8 A2=1.1043D-10 B=5.707D-6 B1=-2.051D-8 C=1.9898D-4 C1=-2.376D-6 D=1.83D-11 E=-0.765D-8 ST=T-273.15 COMP=1.D0-(P/T)*(A+A1*TT+A2*TT**2+(B+B1*TT)*XW * +(C+C1*TT)*XW**2) + (P/T)**2 * (D+EXW**2) RETURN END C C********************************************************************* C