c     program ray1.f computes the rayleigh optical thickness of the standard air
c
      implicit real*8 (a-h,o-z)
c
      real*8 xlamb(2251),xns(2251),zns(2251),xnu(2251),
     1       sig(2251),tr(2251)
      real*8 ylamb(14),yrhon(14),zrhon(2251),rrhon(2251)
c
      data xntot/2.54743d19/
      data a0,a1,a2/6432.8d0,2949810.0d0,25540.0d0/
      data ylamb/0.39d0,0.40d0,0.45d0,0.50d0,0.55d0,0.60d0,0.65d0,
     1           0.70d0,0.75d0,0.80d0,0.85d0,0.90d0,0.95d0,1.0d0/
      data yrhon/1.051d0,1.051d0,1.050d0,1.049d0,1.049d0,1.048d0,
     1           1.048d0,1.048d0,1.048d0,1.047d0,1.047d0,1.047d0,
     2           1.047d0,1.047d0/
c
      nlamb=2251
c     fill the wavelength buffer
      xlamb(1)=0.200 
      do i=1,nlamb-1
         if(i.le.nlamb-1)then
            xlamb(i+1)=xlamb(i)+0.0010
         else
            xlamb(i+1)=0.0d0
         endif
      enddo


      open(unit=8,file='../lores/rayleigh_zia.txt',status='unknown')
      pi=dacos(-1.0d0)
c     determine rhon
      do i=1,nlamb
         if(xlamb(i).le.ylamb(1))then
            zrhon(i)=yrhon(1)
            goto 200
         endif
         if(xlamb(i).gt.ylamb(14))then
            zrhon(i)=yrhon(14)
            goto 200
         endif
         do j=1,13
            if(xlamb(i).gt.ylamb(j) .and. xlamb(i).le.ylamb(j+1))then
              slope=(yrhon(j)-yrhon(j+1))/(ylamb(j)-ylamb(j+1))
              zrhon(i)=yrhon(j)+slope*(xlamb(i)-ylamb(j))
              goto 100
             endif
         enddo
100      continue
200      continue
      enddo
 
      do i=1,nlamb
         xnu(i)=1.0d0/xlamb(i)
         temp=(a0+a1/(146.0d0-xnu(i)**2)+a2/(41.0d0-xnu(i)**2) )
         xns(i)=temp*1.0d-8+1.0d0
         zns(i)=(xns(i)-1.0d0)
         rhon=(6.0d0*(zrhon(i)-1.0d0))/(3.0d0+7.0d0*zrhon(i))
c
c using fixed dpol
c
c        rhon=0.0279d0
c        rhon=0.031d0
c
         rrhon(i)=rhon
         factor=(6.0d0+3.0d0*rhon)/(6.0d0-7.0d0*rhon)
         const=(8.0d0*pi**3)/3.0d0
         bbn=(xns(i)**2-1.0d0)**2
         bbd=(xlamb(i)*1.0d-4)**4*xntot**2
         bb=bbn/bbd
         sig(i)=const*bb*factor
         write(*,*) xlamb(i),rhon,0.0279d0
      enddo
c     assume total number of particles in a verticle column=aaa
c     =2.149d+25 (see penndorf j. opt. soc. am. vol47,176-182,1957)
      aaa=2.149d+25
      do i=1,nlamb
         tr(i)=sig(i)*aaa
      enddo
c
c      write(8,'(a)')'  no          wavl       trr       rhon'
      do i=1,nlamb
c         write(8,160)i-1,xlamb(i),tr(i),rrhon(i)
c160      format(i5,5x,f10.5,1x,f10.5,f10.5)
         write(8,160) xlamb(i)*1000,tr(i)
160      format(f10.1,1x,f10.5)
      enddo

c
c
      stop
      end

