Initial commit

1 files changed, 343 insertions, 0 deletions

diff --git a/synthe/broadenx.for b/synthe/broadenx.for
new file mode 100644
index 0000000..28e057c
--- /dev/null
+++ b/synthe/broadenx.for
@@ -0,0 +1,343 @@
+      PROGRAM BROADENX
+c     constants given D exponents
+c     revised 31jan01 changed read statement for wledge
+c     
+C     TAPE5=INPUT
+C     TAPE6=OUTPUT
+C     TAPE21=SPECTRUM INPUT
+C     TAPE22=SPECTRUM OUTPUT
+C     THE MINIMUM DIMENSION OF H IS (NWL+19999+19999)*NMU
+C     FOR FLUX SPECTRA NMU IS 1
+C      DIMENSION H(4000000)
+      DIMENSION H(4100000)
+      DIMENSION RED(20000),BLUE(20000)
+      DIMENSION RED1(20000),BLUE1(20000),RED2(20000),BLUE2(20000)
+      EQUIVALENCE (RED(1),RED1(1)),(BLUE(1),BLUE1(1))
+      CHARACTER*10 A,B
+      CHARACTER*50 COMMENT
+      REAL*8 LINDAT8(14)
+      REAL*4 LINDAT4(28)
+      DIMENSION XMU(20),QMU(40),WLEDGE(377),TITLE(74)
+      REAL*8 TEFF,GLOG,TITLE,WBEGIN,RESOLU,XMU,WLEDGE,RATIO
+      REAL*8 QMU
+      CHARACTER*1 APLOT(101)
+      DATA APLOT/101*' '/
+      LINOUT=300
+      REWIND 21
+      READ(21)TEFF,GLOG,TITLE,WBEGIN,RESOLU,NWL,IFSURF,NMU,XMU,NEDGE,
+     1(WLEDGE(iedge),iedge=1,nedge)
+      WRITE(6,1010)TEFF,GLOG,TITLE
+ 1010 FORMAT(  5H TEFF,F7.0,7H   GRAV,F7.3/7H TITLE ,74A1)
+      WRITE(6,1007)NMU,(XMU(IMU),IMU=1,NMU)
+ 1007 FORMAT(I4,20F6.3)
+C     IFSURF=3 FOR ROTATED SPECTRUM
+      IF(IFSURF.EQ.3) NMU=1
+      RATIO=1.+1./RESOLU
+      WEND=WBEGIN*RATIO**(NWL-1)
+      WCEN=(WBEGIN+WEND)*.5
+      VSTEP=2.99792458D5/RESOLU
+      WRITE(6,1005)WBEGIN,WEND,RESOLU,VSTEP,NWL
+ 1005 FORMAT(2F14.5,F12.1,F12.5,I10)
+      NMU1=NMU+1
+      NMU2=NMU+NMU
+C
+C     SAMPLE CARDS RIGHT SHIFTED BY 1
+C     THESE ARE EVALUATED AT WLCEN AND HAVE CONSTANT 
+C     RESOLUTION OR RESOLVING POWER
+CGAUSSIAN  3.5       KM                  COMMENT FIELD
+CGAUSSIAN  100000.   RESOLUTION          COMMENT FIELD
+CGAUSSIAN  7.        PM                  COMMENT FIELD
+CGAUSSIAN  .01       CM-1                COMMENT FIELD
+CSINX/X    3.5       KM                  COMMENT FIELD
+CSINX/X    100000.   RESOLUTION          COMMENT FIELD
+CSINX/X    7.        PM                  COMMENT FIELD
+CSINX/X    .01       CM-1                COMMENT FIELD
+CRECT      7.        PM                  COMMENT FIELD
+CRECT      3.5       KM                  COMMENT FIELD
+CRECT      100000.   RESOLUTION          COMMENT FIELD
+CREDT      .01       CM-1                COMMENT FIELD
+CMACRO     2.0       KM                  COMMENT FIELD
+CPROFILE   5.        POINTS              COMMENT FIELD
+CRED       .3        .1        .1        .1        .05
+CBLUE      .3        .1        .1        .1        .05
+C     POINTS ARE TABULATED AT THE SPACING OF THE COMPUTED SPECTRUM
+C     THE CENTER IS THE FIRST POINT FOR EACH WING
+C     TAKING THE CENTER ONLY ONCE, THE PROFILE SUMS TO 1.
+C
+C     IN THESE PROFILES THE RESOLUTION OR RESOLVING POWER VARIES
+C     THEY ARE LINEARLY INTERPOLATED FROM WLBEG TO WLEND
+C          AT WLBEG            AT WLEND
+CGAUSSIAN  3.5       KM        3.        COMMENT FIELD
+CGAUSSIAN  100000.   RESOLUTION 120000.  COMMENT FIELD
+CGAUSSIAN  7.        PM        7.        COMMENT FIELD
+CGAUSSIAN  .01       CM-1      .01       COMMENT FIELD
+CSINX/X    3.5       KM        3.        COMMENT FIELD
+CSINX/X    100000.   RESOLUTION 120000.  COMMENT FIELD
+CSINX/X    7.        PM        7.        COMMENT FIELD
+CSINX/X    .01       CM-1      .01       COMMENT FIELD
+CRECT      7.        PM        7.        COMMENT FIELD
+CRECT      3.5       KM        3.        COMMENT FIELD
+CRECT      100000.   RESOLUTION 120000.  COMMENT FIELD
+CRECT      .01       CM-1      .01       COMMENT FIELD
+CPROFILE   5.        POINTS    5.        COMMENT FIELD
+CRED       .3        .1        .1        .1        .05
+CBLUE      .3        .1        .1        .1        .05
+CRED       .25       .125      .1        .1        .05
+CBLUE      .25       .125      .1        .1        .05
+
+      READ(5,1)A,X1,B,X2,COMMENT
+    1 FORMAT(A10,F10.2,A10,F10.2,A50)
+      WRITE(6,2)A,X1,B,X2,COMMENT
+    2 FORMAT(1X,A10,F10.2,A10,F10.2,A50)
+      FWHM=-1.
+      IF(B.EQ.'PM        ')FWHM=X1/WCEN/1000.*299792.458D0
+      IF(B.EQ.'KM        ')FWHM=X1
+      IF(B.EQ.'RESOLUTION')FWHM=299792.458D0/X1
+      IF(B.EQ.'CM-1      ')FWHM=X1/(1.E7/WCEN)*299792.458D0
+      IF(FWHM.LT.0.)THEN
+      WRITE(6,3)
+      CALL EXIT
+      ENDIF
+      FWHM1=FWHM
+      FWHM2=FWHM
+      IF(X2.GT.0.)THEN
+      IF(B.EQ.'PM        ')FWHM2=X2/WCEN/1000.*299792.458D0
+      IF(B.EQ.'KM        ')FWHM2=X2
+      IF(B.EQ.'RESO      ')FWHM2=299792.458D0/X2
+      IF(B.EQ.'CM-1      ')FWHM2=X2/(1.E7/WCEN)*299792.458D0
+      ENDIF
+      IF(A.EQ.'MACRO     ')GO TO 10
+      IF(A.EQ.'GAUSSIAN  ')GO TO 20
+      IF(A.EQ.'SINX/X    ')GO TO 60
+      IF(A.EQ.'RECT      ')GO TO 30
+      IF(A.EQ.'PROFILE   ')GO TO 40
+      WRITE(6,3)
+    3 FORMAT(10H0BAD INPUT)
+      CALL EXIT
+C
+C     MACROTURBULENT VELOCITY IN KM
+   10 VMAC=X1
+      DO 11 I=1,20000
+      RED(I)=EXP(-(FLOAT(I-1)*VSTEP/VMAC)**2)
+      IF(RED(I).LT.1.E-5)GO TO 12
+   11 CONTINUE
+   12 NPROF=I
+      RED(1)=RED(1)/2.
+      SUM=0.
+      DO 13 I=1,NPROF
+   13 SUM=SUM+RED(I)
+      SUM=SUM*2.
+      DO 14 I=1,NPROF
+      RED(I)=RED(I)/SUM
+   14 BLUE(I)=RED(I)
+      GO TO 50
+C
+C     GAUSSIAN INSTRUMENTAL PROFILE HALF WIDTH IN KM  FWHM
+   20 DO 21 I=1,20000
+      RED(I)=EXP(-(FLOAT(I-1)*VSTEP/FWHM*.8325546D0*2.)**2)
+      IF(RED(I).LT.1.E-5)GO TO 22
+   21 CONTINUE
+   22 NPROF=I
+      RED(1)=RED(1)/2.
+      SUM=0.
+      DO 23 I=1,NPROF
+   23 SUM=SUM+RED(I)
+      SUM=SUM*2.
+      DO 24 I=1,NPROF
+      RED(I)=RED(I)/SUM
+   24 BLUE(I)=RED(I)
+      IF(X2.EQ.0.)GO TO 50
+      DO 25 I=1,20000
+      RED2(I)=EXP(-(FLOAT(I-1)*VSTEP/FWHM2*.8325546D0*2.)**2)
+      IF(RED2(I).LT.1.E-5)GO TO 26
+   25 CONTINUE
+   26 NPROF=I
+      RED2(1)=RED2(1)/2.
+      SUM=0.
+      DO 27 I=1,NPROF
+   27 SUM=SUM+RED2(I)
+      SUM=SUM*2.
+      DO 28 I=1,NPROF
+      RED2(I)=RED2(I)/SUM
+   28 BLUE2(I)=RED2(I)
+      GO TO 50
+C
+C     SINX/X INSTRUMENTAL PROFILE HALF WIDTH IN KM  FWHM
+C     APODIZED BY EXP(-0.06*X**2)
+   60 RED(1)=0.5
+      DO 61 I=2,20000
+      X=(FLOAT(I-1)*VSTEP/FWHM*2.*1.8954942D0)
+      RED(I)=SIN(X)/X*EXP(-0.06*X**2)
+      IF(ABS(RED(I))+ABS(RED(I-1)).LT.1.E-5)GO TO 62
+   61 CONTINUE
+   62 NPROF=I
+      SUM=0.
+      DO 63 I=1,NPROF
+   63 SUM=SUM+RED(I)
+      SUM=SUM*2.
+      DO 64 I=1,NPROF
+      RED(I)=RED(I)/SUM
+   64 BLUE(I)=RED(I)
+      IF(X2.GT.0)GO TO 50
+      RED2(1)=0.5
+      DO 65 I=2,20000
+      X=(FLOAT(I-1)*VSTEP/FWHM2*2.*1.8954942D0)
+      RED2(I)=SIN(X)/X*EXP(-0.06*X**2)
+      IF(ABS(RED2(I))+ABS(RED2(I-1)).LT.1.D-5)GO TO 66
+   65 CONTINUE
+   66 NPROF=I
+      SUM=0.
+      DO 67 I=1,NPROF
+   67 SUM=SUM+RED2(I)
+      SUM=SUM*2.
+      DO 68 I=1,NPROF
+      RED2(I)=RED2(I)/SUM
+   68 BLUE2(I)=RED2(I)
+      GO TO 50
+C
+C     RECTANGULAR INSTRUMENTAL PROFILE HALF WIDTH IN KM  FWHM
+   30 XRECT=FWHM/2./VSTEP
+      NRECT=XRECT+1.5
+      NPROF=NRECT
+      DO 31 I=1,NPROF
+   31 RED(I)=1.
+      RED(NPROF)=XRECT+1.5-FLOAT(NRECT)
+      RED(1)=RED(1)/2.
+      SUM=0.
+      DO 33 I=1,NPROF
+   33 SUM=SUM+RED(I)
+      SUM=SUM*2.
+      DO 34 I=1,NPROF
+      RED(I)=RED(I)/SUM
+   34 BLUE(I)=RED(I)
+      IF(X2.EQ.0.)GO TO 50
+      XRECT=FWHM2/2./VSTEP
+      NRECT=XRECT+1.5
+      NPROF=NRECT
+      DO 35 I=1,NPROF
+   35 RED2(I)=1.
+      RED2(NPROF)=XRECT+1.5-FLOAT(NRECT)
+      RED2(1)=RED2(1)/2.
+      SUM=0.
+      DO 36 I=1,NPROF
+   36 SUM=SUM+RED2(I)
+      SUM=SUM*2.
+      DO 37 I=1,NPROF
+      RED2(I)=RED2(I)/SUM
+   37 BLUE2(I)=RED2(I)
+      GO TO 50
+C
+C     INSTRUMENTAL PROFILE TABULATED AT SPECTRUM POINT SPACING.
+C        RED AND BLUE WINGS BOTH START AT CENTRAL POINT.
+C        THE PROFILE SHOULD SUM TO 1.
+   40 NPROF=X1
+      READ(5,41)(RED(I),I=1,NPROF)
+      READ(5,41)(BLUE(I),I=1,NPROF)
+   41 FORMAT(10X,7F10.6)
+      RED(1)=RED(1)/2.
+      BLUE(1)=BLUE(1)/2.
+      IF(X2.LE.0.)GO TO 50
+      READ(5,41)(RED2(I),I=1,NPROF)
+      READ(5,41)(BLUE2(I),I=1,NPROF)
+      RED2(1)=RED2(1)/2.
+      BLUE2(1)=BLUE2(1)/2.
+   50 WRITE(6,51)(I,RED1(I),BLUE1(I),RED2(I),BLUE2(I),I=1,NPROF)
+   51 FORMAT(I5,4F10.6)
+      WRITE(22)TEFF,GLOG,TITLE,WBEGIN,RESOLU,NWL,IFSURF,NMU,XMU,NEDGE,
+     1WLEDGE
+      NH=(NWL+19999+19999)*NMU
+      DO 52 I=1,NH
+   52 H(I)=0.
+      WRITE(6,1117)
+ 1117 FORMAT(1H1)
+      IF(NMU.EQ.1)GO TO 150
+      DO 57 IWL=1,NWL
+      READ(21)(QMU(I),I=1,NMU)
+      DO 56 IMU=1,NMU
+      Q=QMU(IMU)
+      IWL1000=(IWL+20000)*NMU+IMU
+      DO 53 I=1,NPROF
+   53 H(IWL1000-I*NMU)=H(IWL1000-I*NMU)+BLUE(I)*Q
+      IWL998=(IWL+19998)*NMU+IMU
+      DO 54 I=1,NPROF
+   54 H(IWL998+I*NMU)=H(IWL998+I*NMU)+RED(I)*Q
+   56 CONTINUE
+   57 CONTINUE
+      IF(X2.EQ.0.)GO TO 160
+      REWIND 21
+      READ(21)
+      DO 58 IWL=1,NWL
+      IWLNMU=(IWL+19999)*NMU
+      WT2=FLOAT(IWL-1)/FLOAT(NWL-1)
+      WT1=1.-WT2
+      DO 58 IMU=1,NMU
+   58 H(IWLNMU+IMU)=H(IWLNMU+IMU)*WT1
+      DO 257 IWL=1,NWL
+      READ(21)(QMU(I),I=1,NMU)
+      WT2=FLOAT(IWL-1)/FLOAT(NWL-1)
+      DO 256 IMU=1,NMU
+      Q=QMU(IMU)*WT2
+      IWL1000=(IWL+20000)*NMU+IMU
+      DO 253 I=1,NPROF
+  253 H(IWL1000-I*NMU)=H(IWL1000-I*NMU)+BLUE(I)*Q
+      IWL998=(IWL+19998)*NMU+IMU
+      DO 254 I=1,NPROF
+  254 H(IWL998+I*NMU)=H(IWL998+I*NMU)+RED(I)*Q
+  256 CONTINUE
+  257 CONTINUE
+      GO TO 160
+  150 DO 157 IWL=1,NWL
+      READ(21)QMU(1)
+      Q=QMU(1)
+      IWL1001=IWL+20001
+      DO 153 I=1,NPROF
+  153 H(IWL1001-I)=H(IWL1001-I)+BLUE(I)*Q
+      IWL999=IWL+19999
+      DO 154 I=1,NPROF
+  154 H(IWL999+I)=H(IWL999+I)+RED(I)*Q
+  157 CONTINUE
+      IF(X2.EQ.0.)GO TO 160
+      REWIND 21
+      READ(21)
+      DO 358 IWL=1,NWL
+      WT2=FLOAT(IWL-1)/FLOAT(NWL-1)
+      WT1=1.-WT2
+  358 H(IWL+20000)=H(IWL+20000)*WT1
+      DO 357 IWL=1,NWL
+      READ(21)QMU(1)
+      WT2=FLOAT(IWL-1)/FLOAT(NWL-1)
+      Q=QMU(1)*WT2
+      IWL1001=IWL+20001
+      DO 353 I=1,NPROF
+  353 H(IWL1001-I)=H(IWL1001-I)+BLUE(I)*Q
+      IWL999=IWL+19999
+      DO 354 I=1,NPROF
+  354 H(IWL999+I)=H(IWL999+I)+RED(I)*Q
+  357 CONTINUE
+  160 REWIND 21
+      READ(21)
+      DO 170 IWL=1,NWL
+      READ(21)(QMU(IMU),IMU=1,NMU2)
+      IWLNMU=(IWL+19999)*NMU
+      DO 162 IMU=1,NMU
+  162 QMU(IMU)=H(IWLNMU+IMU)
+      WRITE(22)(QMU(I),I=1,NMU2)
+      IF(IWL.GT.LINOUT)GO TO 170
+      WAVE=WBEGIN*RATIO**(IWL-1)
+      RESID=QMU(1)/QMU(NMU1)
+      IRESID=RESID*1000.+.5
+      IPLOT=RESID*100.+1.5
+      IPLOT=MAX0(1,MIN0(101,IPLOT))
+      APLOT(IPLOT)='X'
+      WRITE(6,2300)IWL,WAVE,IRESID,APLOT
+ 2300 FORMAT(1H ,I5,F11.4,I7,101A1)
+      APLOT(IPLOT)=' '
+  170 CONTINUE
+      READ(21)NLINES
+      WRITE(22)NLINES
+      DO 200 I=1,NLINES
+      READ(21)LINDAT8,LINDAT4
+      WRITE(22)LINDAT8,LINDAT4
+  200 CONTINUE
+      CALL EXIT
+      END