      program spotline
      implicit double precision (a-z)
      character *12 nameout
      integer wt,wt1,wt2,wt3,wt4,wt5,wla,vis,wyd,i,j,k,kk,nul
c     ****************************************************************
      write(*,*) ' program   S P O T L I N E   ver 5.00 (25 XI 1996)'
      write(*,*) '                             by OWal'
      write(*,*) ' start!!'
      write(*,*)
c     CONST. DEFINITION **********************************************
      pi=acos(-1)
      k=64
      kk=65
      pi2=pi/2
      rad=pi/180
      wt1=0
      wt2=0
      wt3=0
      wt4=0
      wt5=0
      wt=1
      wla=0
      vis=1

c     ****************************************************************
c     wt? =1 when tau? exist , and 0 when tau? not exist
c     when spot is visible vis=1, spot not visible vis=0
c     when wla=1 Ilum=Ilum(la cetral.)
c     dok - level of machine error
c     nul 0 - all tau   1 - rezults only if spot is visible
c     READING DATA ***************************************************
      open(1,file='spotline.cfg',status='old',err=997)
      read(1,*,err=998) al
      read(1,*,err=998) be
      read(1,*,err=998) inkl
      read(1,*,err=998) lao
      read(1,*,err=998) lastart
      read(1,*,err=998) lastop
      read(1,*,err=998) ila
      read(1,*,err=998) vo
      read(1,*,err=998) u
      read(1,*,err=998) c
      read(1,*,err=998) itau
      read(1,*,err=998) dok
      read(1,*,err=998) nul
      read(1,*,err=998) wyd
      read(1,*,err=995) tau
123   write(1,887)
      write(1,*) '   tau,   tau1,   tau2,   tau3,   tau4,   tau5,   tau6
     #,   la1min,    la2min,    la2max,  la1max,   DeltaLa, vis,   file,
     #   La.cen,     I.cen'
      write(1,887)

      stau=1/itau
      sla=(lastop-lastart)/ila
      al=al*rad
      be=be*rad
      inkl=inkl*rad
      a=(c/vo)

c     CHECKING A CONDITIOONS FOR TAU *********************************
      if((be-inkl-al).ge.pi2) goto 996

c     TAU 1  ***********************
      if(((al+abs(be-inkl)).lt.pi2).and.(pi2.lt.(al+be+inkl))) then
          wt1=1
          tau1=acos((sin(al)-cos(be)*cos(inkl))/(sin(be)*sin(inkl)))
          tau1=tau1/(2*pi)
          write(*,*) 'tau1=',tau1
        else
          wt1=0
      end if

C     TAU 2  ***********************
      if(((be-inkl-al).lt.pi2).and.(pi2.lt.(be+inkl-al))) then
          tau2=acos((-1*sin(al)-cos(be)*cos(inkl))/(sin(be)*sin(inkl)))
          tau2=tau2/(2*pi)
          if(wt1.eq.0)  tau1=0
          wt2=1
          write(*,*) 'tau2=',tau2
        else
          if(wt1.eq.1)  tau2=0.6
          wt2=0
      end if

C     TAU 3  ***********************
      if((pi2.lt.(be+inkl)).and.((be+inkl).lt.pi).and.
     # (abs(be-inkl).lt.pi2)) then
          tau3=acos(-1/tan(inkl)/tan(be))
          tau3=tau3/(2*pi)
          wt3=1
          write(*,*) 'tau3=',tau3
        else
          wt3=0
      end if

C     TAU 4  ***********************
      if(pi2.lt.(al+be)) then
      pom=cos(al)/sin(be)
      if(pom.eq.1) pom=pom-dok
          tau4=asin(pom)
          tau4=tau4/(2*pi)
          wt4=1
          tau6=0.5-tau4
          write(*,*) 'tau4=',tau4
          write(*,*) 'tau6=',tau6
        else
          wt4=0
      end if

C     TAU 5  ***********************
      if(inkl.lt.be) then
          tau5=acos(tan(inkl)/tan(be))
          tau5=tau5/(2*pi)
          wt5=1
          write(*,*) 'tau5=',tau5
        else
          wt5=0
      end if

c     POLAR POINT IS IN SPOT AREA ************************************
      if(be.ge.al) then
        if(((be+inkl-al).lt.pi2).and.(pi2.lt.(al+be-inkl))) then
        tau1=0
        tau2=0.6
        endif

        if(((be+inkl-al).lt.pi2).and.(pi2.lt.(inkl-be+al))) then
        tau1=0
        tau2=0.6
        endif
      end if

17    if((al+be+inkl).le.pi2) then
          if((wt1.eq.0).and.(wt2.eq.0)) then
            tau1=0.6
            tau2=0.6
          end if
        else
          if((wt1.eq.0).and.(wt2.eq.0)) then
            tau1=0
            tau2=0.6
          end if
      endif

      if(wt3.eq.0) tau3=tau2
      if(wt4.eq.0) tau4=tau3
      if(wt4.eq.0) tau6=tau2
      if(wt5.eq.0) tau5=tau4

c     TAU LOOP BEGIN *************************************************
50    continue
      do 100 i=0,itau,1
      if(wt.eq.1) goto 124
      tau=-0.5+i*stau
124   la1min=0
      la1max=0
      la2min=0
      la2max=0
      lacen=0
      wla=0
      z1=0
      z2=0
      z3=0
      z4=0

      if(abs(tau).lt.dok) tau=0
      write(*,*) 'tau=',real(tau)

      k=k+1
      if(k.ge.90) then
      kk=kk+1
      k=65
      endif
      vis=1

c     spot is invisible in this tau  *******************************
      if(tau2.lt.abs(tau)) then
        vis=0
        write(*,*) real(tau),'spot is not visible'
        ilum=0
        k=k-1
        nameout='no-file'
        if(nul.ne.0) goto 29
        k=k+1
      endif

      nameout=char(kk)//char(k)//'_out.dat'
      open(2,file=nameout)

29    if(vis.eq.0) goto 30

c     WHEN SPOT IS VISIBLE THEN: *************************************
      vis=1
      chi=acos(sin(be)*sin(tau*2*pi))
      phi=acos(cos(be)*sin(inkl)-sin(be)*cos(inkl)*cos(tau*2*pi))
      pomoc=cos(be)*cos(inkl)+sin(be)*sin(inkl)*cos(tau*2*pi)
      if(pomoc.ge.1.0) pomoc=pomoc-dok
      psi=acos(pomoc)

      la1max=lao*(1-(sin(inkl)*cos(al+chi))/a)
      la1min=lao*(1-(sin(inkl)*cos(al-chi))/a)

c     IF ONLY A PART OF SPOT IS VISIBLE ******************************
      if(abs(tau).lt.tau1) goto 20
      pomoc=sin(al)*sin(al)-cos(psi)*cos(psi)
      if(abs(pomoc).lt.(dok)) pomoc=0

      la2max=cos(al)*cos(chi)-cos(phi)*sqrt(pomoc)
      la2max=sin(inkl)*la2max/(a*sin(psi)*sin(psi))
      la2max=lao*(1-la2max)

      la2min=cos(al)*cos(chi)+cos(phi)*sqrt(pomoc)
      la2min=sin(inkl)*la2min/(a*sin(psi)*sin(psi))
      la2min=lao*(1-la2min)
      goto 30

c     WHOLE SPOT IS VISIBLE ******************************************
20    continue
      la2min=0
      la2max=0

30    continue
      dela=real(abs(la1max-la1min))


c     LAMBDA LOOP BEGIN **********************************************
2000  continue

            x1min=a*(1-(la1min/lao))/sin(inkl)
            x1max=a*(1-(la1max/lao))/sin(inkl)
            x2min=a*(1-(la2min/lao))/sin(inkl)
            x2max=a*(1-(la2max/lao))/sin(inkl)
      if(abs(x1min).gt.abs(x1max)) then
      pom=x1min
      x1min=x1max
      x1max=pom
      end if
      if(abs(x2min).gt.abs(x2max)) then
      pom=x2min
      x2min=x2max
      x2max=pom
      end if
      pom=0
      if(abs(x1min).eq.abs(x1max)) then
      x1min=-1*abs(x1min)
      x1max=abs(x1max)
      end if
      if(abs(x2min).eq.abs(x2max)) then
      x2min=-1*abs(x2min)
      x2max=abs(x2max)
      end if

      do 200 j=0,ila,1
      Icen=0
      la=lastart+j*sla
      if(wla.eq.1) la=lacen
      x=a*(1-(la/lao))/sin(inkl)
      if(vis.eq.0) then
        ilum=0
          z1=0
          z2=0
        goto 40
      endif

c     IF LAMBDA<la1min  or la1max<lambda  I(lambda)=0 ****************
c      lambda poza plama <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      if((la1min.lt.la).and.(la.lt.la1max)) goto 35
      ilum=0
      vis=0
      goto 39

35    continue
      vis=1
      if((wt5.eq.0).and.(wt4.eq.0)) tau5=tau1
      if((wt5.eq.0).and.(wt4.eq.1)) tau5=tau4
      if(wt3.eq.0) tau3=tau2

c     for all lambda *****************
      if((0.le.abs(tau)).and.(abs(tau).lt.abs(tau1))) then
         z1=15
         z2=16
         goto 1000
      end if

c     IF  la1min < lambda < la2min ****************
c     wewn   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

      if( ((x1min.lt.x).and.(x.lt.x2min)).or.
     #((x1min.gt.x).and.(x.gt.x2min))) then
      z1=15
      z2=16
        if((abs(tau3).lt.abs(tau)).and.(abs(tau).lt.abs(tau2))) then
           if(wt3.eq.0) goto 1002
           ilum=0
           z1=0
           z2=0
           goto 40
        end if

        if((abs(tau6).lt.abs(tau)).and.(abs(tau).lt.abs(tau2))) then
           if(wt4.eq.0) goto 1002
           ilum=0
           z1=0
           z2=0
           goto 40
        end if
      endif
1002  continue


c     IF  la2min < lambda < la2max ****************
c     srodek            <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      if(((x2min.lt.x).and.(x.lt.x2max)).or.
     #((x2min.gt.x).and.(x.gt.x2max))) then
      z1=15
      z2=17

        if (pi2.lt.(al+be)) then
           z1=18
           z2=16
           if((abs(tau5).lt.abs(tau)).and.(abs(tau).lt.abs(tau3))) then
              if(wt5.eq.0) goto 1003
              z1=15
              z2=17
           end if
1003       continue
           if((abs(tau3).lt.abs(tau)).and.(abs(tau).lt.abs(tau2))) then
              if(wt3.eq.0) goto 1101
              z1=15
              z2=17
           end if
1101       continue
           if((abs(tau6).lt.abs(tau)).and.(abs(tau).lt.abs(tau2))) then
              if(wt4.eq.0) goto 1102
              z1=15
              z2=17
           end if
1102       continue

           if((abs(tau4).lt.abs(tau)).and.(abs(tau).lt.abs(tau6))) then
              if(wt4.eq.0) goto 1112
              z1=15
              z2=17
           end if
1112       continue

           if((abs(tau4).lt.abs(tau)).and.(abs(tau).lt.abs(tau5))) then
              if(wt4.eq.0) goto 1119
              z1=18
              z2=16
           end if
1119       continue


           if((abs(tau1).lt.abs(tau)).and.(abs(tau).lt.abs(tau4))) then
              if(wt4.eq.0) goto 1113
              z1=15
              z2=17
           end if
1113       continue

        end if
        goto 1000
      end if

c     IF  la2max < lambda < la1max ****************
c      zewn. <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      if(((x2max.lt.x).and.(x.lt.x1max)).or.
     #((x2max.gt.x).and.(x.gt.x1max))) then
        z1=15
        z2=16
        if((abs(tau3).lt.abs(tau)).and.(abs(tau).lt.abs(tau2))) then
          if(wt3.eq.0) goto 1005
          ilum=0
          z1=0
          z2=0
          goto 40
        end if
1005    continue
        if (pi2.lt.(al+be)) then
           z1=18
           z2=17
           if((abs(tau1).le.abs(tau)).and.(abs(tau).le.abs(tau4))) then
              if(wt4.eq.0) goto 1001
              z1=15
              z2=16
           end if
1001       continue
           if((abs(tau3).le.abs(tau)).and.(abs(tau).le.abs(tau2))) then
             if(wt3.eq.0) goto 1104
             ilum=0
             z1=0
             z2=0
             goto 40
           end if
1104       continue
           if((abs(tau6).le.abs(tau)).and.(abs(tau).le.abs(tau2))) then
             if(wt4.eq.0) goto 1000
             ilum=0
             z1=0
             z2=0
             goto 40
           end if

           if((abs(tau1).le.abs(tau)).and.(abs(tau).le.abs(tau6))) then
              if(wt4.eq.0) goto 1006
              z1=15
              z2=16
           end if
1006       continue
           if((abs(tau4).le.abs(tau)).and.(abs(tau).le.abs(tau6))) then
              if(wt4.eq.0) goto 1007
              z1=18
              z2=17
           end if
1007       continue

        end if
      end if

1000  continue
      if(z1.eq.15) then
      z1=sin(al)*sin(al)*sin(chi)*sin(chi)
      z1=z1-(cos(al)*cos(chi)-x)*(cos(al)*cos(chi)-x)
      if(abs(z1).lt.dok) z1=0
      z1=sqrt(z1)
      z1=(cos(phi)*(cos(al)-x*cos(chi))-cos(psi)*z1)
     #/(sin(chi)*sin(chi))
      end if

      if(z2.eq.16) then
      z2=sin(al)*sin(al)*sin(chi)*sin(chi)
      z2=z2-(cos(al)*cos(chi)-x)*(cos(al)*cos(chi)-x)
      if(abs(z2).lt.dok) z2=0
      z2=sqrt(z2)
      z2=(cos(phi)*(cos(al)-x*cos(chi))+cos(psi)*z2)
     #/(sin(chi)*sin(chi))
      endif

      if(z1.eq.17) then
      z1=1-x*x
      if(abs(z1).lt.dok) z1=0
      z1=sqrt(z1)
      end if

      if(z2.eq.17) then
      z2=1-x*x
      if(abs(z2).lt.dok) z2=0
      z2=sqrt(z2)
      end if

      if(z1.eq.18) then
      z1=1-x*x
      if(abs(z1).lt.dok) z1=0
      z1=-1*sqrt(z1)
      end if

      if(z2.eq.18) then
      z2=1-x*x
      if(abs(z2).lt.dok) z2=0
      z2=-1*sqrt(z2)
      end if

      z3=1-x*x
      if(abs(z3).lt.dok) z3=0
      z3=sqrt(z3)
      z4=-1*sqrt(z3)


c     ilumination= ***************************************************
      pom=1-x*x
      if(abs(pom).lt.dok) pom=0
      ilum=(z2/sqrt(pom))
      if(ilum.eq.1) ilum=1-dok
      if(ilum.eq.-1) ilum=1+dok
      ilum=asin(ilum)
      pomoc=(z1/sqrt(pom))
      if(pomoc.eq.1) pomoc=1-dok
      ilum=ilum-asin(pomoc)
      ilum=(1-u)*(z2-z1)+u*(pom)*ilum/2
      pomoc=(1-(x*x)-(z2*z2))
      if(abs(pomoc).lt.dok) pomoc=0
      pom=(1-(x*x)-(z1*z1))
      if(abs(pom).lt.dok) pom=0
      ilum=ilum+u*(z2*sqrt(pomoc)-z1*sqrt(pom))/2
      if(wla.eq.1) icen=ilum
39    vis=1
40    continue

      if(wla.eq.1) goto 3000
      if(abs(x).lt.1) then
      if((wyd.eq.1).and.(ilum.ne.0)) write(2,8891,err=80) x,la,
     #real(z1),real(z2),real(z3),real(z4)
      if(wyd.eq.2) write(2,8892,err=80) x,la,real(ilum),real(ilum+tau)
     #,real(tau)
      if(wyd.eq.0) write(2,889,err=80) x,la,real(ilum),real(ilum+tau),
     #real(tau),real(z1),real(z2),real(z3),real(z4)
      end if
c     ****************************************************************
      z1=0
      z2=0
      z3=0
      z4=0
      pomoc=0
      pom=0
      ilum=0

80    continue
c     END of lambda loop *********************************************
      z1=0
      z2=0
      z3=0
      z4=0
200   continue
      close (2)
3000  continue

      if(wla.eq.0) then
         lacen=(la1min+la1max)/2
         wla=1
         goto 2000
      end if
      wla=0

      write(1,888) real(tau),real(tau1),real(tau2),real(tau3),
     #real(tau4),real(tau5),real(tau6),real(la1min),
     #real(la2min),real(la2max),real(la1max),dela,vis,nameout,
     #real(lacen),real(icen)


c     END OF TAU LOOP *****************************************
      if(wt.eq.1) goto 125
100   continue
125   write(*,*)  '                                T H E   E N D '
      write(1,887)
c     *********************************************************

887   format(144(1H#))
888   format(1x,f6.2,2x,f6.4,2x,f6.4,2x,f6.4,2x,f6.4,2x,f6.4,
     #2x,f6.4,2x,f9.4,1x,f9.4,1x,f9.4,1x,f9.4,1x,f9.4,1x,i2,
     #1x,a13,1x,f9.4,1x,f9.4)
889   format(1x,f7.4,1x,f9.4,1x,f9.6,1x,f9.6,1x,f7.4,1x,f7.4,
     #1x,f7.4,1x,f7.4,1x,f7.4)
8891  format(1x,f7.4,1x,f9.4,1x,f7.4,1x,f7.4,1x,f7.4,1x,f7.4)
8892  format(1x,f7.4,1x,f9.4,1x,f9.6,1x,f9.6,1x,f7.4)

c     errors information *********************************************
      goto 999
995   wt=0
      goto 123
996   write(*,*) 'spot ALWAYS is NOT visible !!!'
      write(1,*) 'spot ALWAYS is NOT visible !!!'
      goto 999
997   write(*,*) 'I can not open a file "spotline.cfg"'
      goto 999
998   write(*,*) 'bad reading data from file "spotline.cfg"'
      goto 999
999   continue
      end
