CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890

      Program calibracja

      real mass, zenith, dsol, lat
      real tu, xlon, xlat, asol, phi0, V0(10), VV0
      real tau, x, a, b, s1, s2, s3, s4, dmin
      real m(100), lnV(10,100), hour(100), Vd(10,100),V(10,100)
      real h(100), min(100), lambda(10), err(10)
      integer  month, day, year, iday, imon
      integer i, j, nrec, ndetect
      Write(*,*) '--------------------------------------------------'
      Write(*,*) '              PROGRAM CALIB1.1                    '
      Write(*,*) '          Autor Krzysztof Markowicz                '
      Write(*,*) '--------------------------------------------------'
      Write(*,*)
      Write(*,*) 'Program wykonuje kalibracje Fotometru metoda Langleya'
      Write(*,*) 'Program czyta dane z pliku calib.dat '
      Write(*,*) 'Informacje w pliku Readme.txt'
      Write(*,*)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC Read data from text file
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      open(20,file='calib.dat',status='old')
      read(20,*) ndetect, nrec
      read(20,*) (h(i),min(i),(Vd(j,i),V(j,i),j=1,ndetect),i=1,nrec)
      print*, 'Plik kalibracyjny wczytany pomyslnie'
      close(20)
       do i=1,nrec
        hour(i)=h(i)+min(i)/60
        do j=1,ndetect
         lnV(j,i)=alog(V(j,i)-Vd(j,i))
        end do
       end do
      Write(*,*) ' Podaj szerokosc geograficzna stacji'
      Read(*,*) xlat
      Write(*,*) ' Podaj dlugosc geograficzna geograficzna stacji'
      Read(*,*) xlon
      Write(*,*) ' Podaj rok w czasie wykonywania kalibracji'
      Read(*,*) year
      Write(*,*) ' Podaj miesiac'
      Read(*,*) month
      Write(*,*) ' Podaj dzien'
      Read(*,*) day
      Write(*,*) ' Podaj cisnienia na poziomie stacji w hpa'
      Read(*,*) p
      Write(*,*) ' Podaj dlugosc fali Fotometru w nano-metrach'
      do j=1,ndetect
        Write(*,*) 'detektor ',j,' :'
        Read(*,*) lambda(j)
      end do
      Write(*,*)
      Write(*,*) '--------------------------------------------------'
      call varsol(day,month,dsol)

      do i=1,nrec
       tu=hour(i)
       call possol (month,day,tu,xlon,xlat,asol,phi0)
       zenith=asol
       call airmass(zenith, mass)
       m(i)=mass
      end do
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C start interpolation
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      do j=1,ndetect
       s1=0
       s2=0
       s3=0
       s4=0
        do i=1,nrec
         s1=s1+lnV(j,i)*m(i)
         s2=s2+lnV(j,i)
         s3=s3+m(i)
         s4=s4+m(i)*m(i)
        end do
        a=(-nrec*s1+s2*s3)/(s3*s3-nrec*s4)
        b=(s2-a*s3)/nrec
        VV0=b-alog(dsol)
        V0(j)=exp(VV0)
       end do
       do j=1,ndetect
        dmin=0.0
        do i=1,nrec
          dmin=dmin+(m(i)*a+b-lnV(j,i))**2
        end do
        err(j)=dmin/nrec
       end do
       Do j=1,ndetect
        print*, 'Detektor',j
        Write(*,'(" Wspolczynnik kalibracyjny wynosi ",f5.3)'),V0(j)
        write(6,'(" Blad aproksymacji wynosi ",f7.5)'),err(j)
        print *,'------------'
       end do

       open (30,file='lang.dat',status='unknown')
       write(30,*) 'Plik kalibracyjny Fotometru metoda Langleya'
       write(30,'(" Szerokosc geograficzna: ",f5.2)'),xlat
       write(30,'(" Dlugosc geograficzna: ",f5.2)'),xlon
       write(30,*) 'Data kalibracji: ',year, month, day
       write(30,*) 'Liczba detektorow: ', ndetect
       write(30,*) 'Dlugosc fali w nm: ',(lambda(j),j=1,ndetect)
       write(30,*) 'Wspolczynnik kalibracyjny: ',(V0(j),j=1,
     $ ndetect)
       write(30,*) 'Blad aproksymacji: ',(err(j),j=1,ndetect)
       write(30,*), ndetect
       write(30,*), (lambda(j),j=1,ndetect)
       write(30,*), (V0(j),j=1,ndetect)
      close(30)
      end


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC Get airmass from solar zenith angle
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      Subroutine airmass(zenith,mass)
       real a, b, c, mass, z
       a = 0.50572
       b = 6.07995
       c = 1.6364
       el=90-zenith
       z=3.14159*el/180
       mass = 1.0 /(sin(z) + a * (el+ b)**-c )
      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine varsol (jday,month,dsol)
      real dsol,pi,om
      integer jday,month,j end

      if (month.le.2) goto 1
      if (month.gt.8) goto 2
      j=31*(month-1)-((month-1)/2)-2+jday
      goto 3
    1 j=31*(month-1)+jday
      goto 3
    2 j=31*(month-1)-((month-2)/2)-2+jday
    3 pi=2.*acos (0.)
      om=(.9856*float(j-4))*pi/180.
      dsol=1./((1.-.01673*cos(om))**2)
      return
      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC Subroutine from 6S package to get Sun position
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      subroutine possol (month,jday,tu,xlon,xlat,
     a asol,phi0)

      real tu,xlon,xlat,asol,phi0
      integer month,jday,ia,nojour

c solar position (zenithal angle asol,azimuthal angle phi0
c in degrees)
c jday is the number of the day in the month

      ia = 0
      call day_number(jday,month,ia,nojour)
      call pos_fft (nojour, tu, xlon, xlat, asol, phi0)
      if(asol.gt.90) call print_error(
     s 'The sun is not raised')
      return
      end

      subroutine day_number(jday,month,ia,j)
      integer jday, month, ia, j

      if (month.le.2) then
                      j=31*(month-1)+jday
                      return
                      endif
      if (month.gt.8) then
                      j=31*(month-1)-((month-2)/2)-2+jday
                      else
                      j=31*(month-1)-((month-1)/2)-2+jday
                      endif
      if(ia.ne.0 .and. mod(ia,4).eq.0) j=j+1
      return
      end

      subroutine pos_fft (j,tu,xlon,xlat,asol,phi0)
      real tu, xlat, asol,phi0, tsm, xlon,xla, xj, tet,
     a a1, a2, a3, a4, a5, et, tsv, ah, b1, b2, b3, b4,
     a b5, b6, b7, delta, amuzero, elev, az, caz, azim, pi2
      integer j
      parameter (pi=3.14159265,fac=pi/180.)
c solar position (zenithal angleasol,azimuthal angle phi0
c in degrees)
c j is the day number in the year
c
c mean solar time (heure decimale)

      tsm=tu+xlon/15.
      xla=xlat*fac
      xj=float(j)
      tet=2.*pi*xj/365.

c time equation (in mn.dec)
      a1=.000075
      a2=.001868
      a3=.032077
      a4=.014615
      a5=.040849
      et=a1+a2*cos(tet)-a3*sin(tet)-a4*cos(2.*tet)-a5*sin(2.*tet)
      et=et*12.*60./pi

c true solar time

      tsv=tsm+et/60.
      tsv=(tsv-12.)

c hour angle
      ah=tsv*15.*fac

c solar declination (in radian)
      b1=.006918
      b2=.399912
      b3=.070257
      b4=.006758
      b5=.000907
      b6=.002697
      b7=.001480
      delta=b1-b2*cos(tet)+b3*sin(tet)-b4*cos(2.*tet)+b5*sin(2.*tet)-
     &b6*cos(3.*tet)+b7*sin(3.*tet)

c elevation,azimuth

      amuzero=sin(xla)*sin(delta)+cos(xla)*cos(delta)*cos(ah)
      elev=asin(amuzero)
      az=cos(delta)*sin(ah)/cos(elev)
      if ( (abs(az)-1.000).gt.0.00000) az = sign(1.,az)
      caz=(-cos(xla)*sin(delta)+sin(xla)*cos(delta)*cos(ah))/cos(elev)
      azim=asin(az)
      if(caz.le.0.) azim=pi-azim
      if(caz.gt.0.and.az.le.0) azim=2*pi+azim
      azim=azim+pi
      pi2=2*pi
      if(azim.gt.pi2) azim=azim-pi2
      elev=elev*180./pi
c conversion in degrees
      asol=90.-elev
      phi0=azim/fac
      return
      end

      subroutine print_error(tex)
      character *(*) tex
      logical ier
      integer iwr
c JM 6/99 common/sixs_ier/iwr,ier
      iwr = 6
      ier = .TRUE.
      write(iwr,'(a)')tex
      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC