C===========================================================================
C
C	program: aspot.for
C
C	purpose: calculates light curve for star with circular spots
C	see budding, e., 1977.  astr sp. sci., 48, 207.
C
C	subroutines: spot
C
C	input: tty
C
C	output: datin.dat, spot.dat
C
C===========================================================================
C
	implicit real*8 (a-h, o-z)
	dimension wk(10),gam0(10),bet0(10),djmin(10),per(10)
	dimension vbar(10),u(10),v(10)
	character yn*1,col*8,titl*20
C
	write(*,'(
     +''This program computes a theoretical starspot model
     + from input data.''/''Random noise can be included if desired.''/
     +''Output data is written to DATIN.DAT and SPOT.DAT.''/)')
C
	jer = 69093
	write(*,'(''Number of wavelengths? '',$)')
	read(*,*) nwav
C
	do i=1,nwav
	 write(*,'(/''Wavelenth '',i4)') i
	 write(*,'(''Unspotted magnitude? '',$)')
	 read(*,*) vbar(i)
	 write(*,'(''Limb-darkening coefficient? '',$)')
	 read(*,*) u(i)
	enddo
C
	write(*,'(''Fiducial intensity level of photosphere (0-1)?
     + '',$)')
	read(*,*) uf
	write(*,'(''Inclination of rotational axis (deg)? '',$)')
	read(*,*) ync0
	frad = 0.017453293
	ync = frad*ync0
C
	write(*,'(''No of spots? '',$)')
	read(*,*) jum
C
	do i=1,jum
	 write(*,'(''Frequency (-ve for period)? '',$)')
	 read(*,*) fr
	 if(fr.lt.0.0) then
	  per(i) = -fr
	 else
	  per(i) = 1.0/fr
	 endif
	 write(*,'(''HJD of light minimum? '',$)')
	 read(*,*) djmin(i)
	 write(*,'(''Fractional brightness of spot? '',$)')
	 read(*,*) wk(i)
	 write(*,'(''Radius of spot in degrees? '',$)')
	 read(*,*) gam0(i)
	 write(*,'(''Latitude of spot in degrees? '',$)')
	 read(*,*) bet0(i)
	enddo
C
	nv = 0
C
	write(*,'(//''Read HJDs from DATIN.DAT (y/n)? '',$)')
	read(*,1000) yn
1000	format(a)
	if(yn.eq.'y'.or.yn.eq.'Y') then
	 nv = 1
	 call openold( 2)
	 call opennew( 1)
	 read(2,1000) titl
	 write(1,1000) titl
	 read(2,*) p,epo,mv
	 do i=1,mv
 	  read(2,1000) col
	 enddo
	 read(2,1000) titl
	 write(1,1002) p,epo,nwav
1002	 format(f12.6/f14.3/i4)
	 write(1,1003) (u(i),i=1,nwav)
1003	 format(f8.4)
	 write(1,1000) titl
	else
	 call opennew( 1)
	 write(*,'(''Grand epoch? '',$)')
	 read(*,*) epo
	 write(1,1004) epo,nwav
1004	 format('Starspot model'/' 1.0000000'/f12.3/i4)
	 write(1,1003) (u(i),i=1,nwav)
	 write(1,'(''Starspot model'')')
	 write(*,'(''Starting HJD? '',$)')
	 read(*,*) dj0
	 write(*,'(''Stopping HJD? '',$)')
	 read(*,*) dj1
	 write(*,'(''Step in days? '',$)')
	 read(*,*) deljd
	endif
C
	write(*,'(''Random se? '',$)')
	read(*,*) rse
C
C	loop through no of data points
C
	if(nv.eq.0) dj = dj0 - deljd
	n = 0
10	n = n + 1
	if(nv.eq.0) then
	 dj = dj + deljd
	else
12  	 read(2,*,end=20) dj
	 if(dj.lt.0.0) then
	  vv = 0.0
	  write(1,1010) dj,n,(vv,iw=1,nwav)
	  read(2,1000,end=20) titl
	  write(1,1000) titl
	  goto 12
	 endif
	endif
C
C	loop through wavelength
	do iw=1,nwav
	 yy = 0.0
C
C	loop through number of spots
	 do j = 1,jum
	  phi = (dj - djmin(j))/per(j)
	  k = phi
	  phi = phi - k
	  ph = 6.28318*phi
	  gam = frad*gam0(j)
	  bet = frad*bet0(j)
	  call spot(wk(j),u(iw),gam,bet,ync,uf,ph,y)
	  yy = yy + y
	 end do
	 v(iw) = -2.5*log10(1.0-yy) + vbar(iw)
	 if(rse.ne.0.0) call erran(v(iw),rse,jer)
	enddo
C
 	write(1,1010) dj,n,(v(iw),iw=1,nwav)
1010	format(f14.4,i8,10f10.4)
	if(nv.eq.0.and.dj.ge.dj1) goto 20
	goto 10
C
20	close(2)
	open(unit=2,file='spot.dat',status='unknown')
	do i = 1,jum
	 write(2,1020) wk(i),u,gam0(i),bet0(i),ync0,r,vbar,djmin(i)
1020	 format(e16.7)
	end do
C
	close(1)
	close(2)
	call progrun('spot')
	stop
	end
C===========================================================================
C
C		spot(wk,u,gam,bet,ync,r,phi,y)
C
C===========================================================================
C
	subroutine spot(wk,u,gam,bet,ync,r,phi,y)
	implicit real*8 (a-h, o-z)
C
	fac = 3.0*r*(1.0-wk)/(3.0-u)
	sg = sin(gam)
	sg2 = sg*sg
	cg = cos(gam)
	cg2 = cg*cg
	tg2 = sg2/cg2
	t1 = cos(bet)*sin(ync)
	t2 = sin(bet)*cos(ync)
	z0 = t1*cos(phi)+t2
	th0 = acos(z0)
	q1 = cos(th0+gam)
	q2 = cos(th0-gam)
	z1 = sqrt(1.0-z0*z0)
	d = cg*z1
C
C	annular case:
	if(q1.ge.0.0.and.q2.ge.0.0) then
	 s0=sg2*z0
	 s1=0.6666667*(1.0-cg*(cg2+1.50*d*d*tg2))
	end if
C
C	partial case:
 	if(q1.lt.0.0.and.q2.ge.0.0) then
	 s = cg2/d
	 v = (d-s)/(z0*sg)
	 if(abs(s).gt.1.0) then
 	  if(s.lt.0.0) s=-1.0
	  if(s.gt.0.0) s=1.0
	 end if
	 if(abs(v).gt.1.0) then
	  if(v.lt.0.0) v=-1.0
	  if(v.gt.0.0) v=1.0
	 end if
	 vs=v/s
	 if(abs(vs).gt.1.0) then
	  if(vs.lt.0.0) vs=-1.0
	  if(vs.gt.0.0) vs=1.0
	 end if
	 f1 = acos(s)-s*sqrt(1.0-s*s)
	 f2 = acos(v)-v*sqrt(1.0-v*v)
	 s0 = 0.318309886*(f1+sg2*z0*f2)
	 f3 = sg*z0*(3.0*sg2-1.0)*sqrt(1.0-v*v)
 	 f4 = 2.0*s*cg2+3.0*d*sg2
	 f5 = cg/(2.0*s)
	 s1 = (acos(vs)+f5*(f3-f4*acos(v)))*0.212206591
	end if
C
	if(q1.lt.0.0.and.q2.lt.0.0) then
	 s0 = 0.0
	 s1 = 0.0
	end if
C
	y = fac*((1.0-u)*s0+u*s1)
	return
	end
C===========================================================================
C
C		erran(v, se, ij)
C
C	introduces gaussian noise
C
C===========================================================================
C
	subroutine erran(v,se,ij)
	implicit real*8 (a-h, o-z)
C
	s = 0.0
	do i=1,12
	 s = s + ran3(ij)
	end do
 	v = v + se*(s - 6.0)
	return
	end
C
C===============================================================================
C
C				RAN3
C
C	Random number generator
C
C===============================================================================
C
      FUNCTION ran3(idum)
	implicit real*8 (a-h,o-z)
      INTEGER idum
      INTEGER MBIG,MSEED,MZ
      PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG)
      INTEGER i,iff,ii,inext,inextp,k
      INTEGER mj,mk,ma(55)
      SAVE iff,inext,inextp,ma
      DATA iff /0/
C
      if(idum.lt.0.or.iff.eq.0)then
        iff=1
        mj=MSEED-iabs(idum)
        mj=mod(mj,MBIG)
        ma(55)=mj
        mk=1
        do 11 i=1,54
          ii=mod(21*i,55)
          ma(ii)=mk
          mk=mj-mk
          if(mk.lt.MZ)mk=mk+MBIG
          mj=ma(ii)
11      continue
        do 13 k=1,4
          do 12 i=1,55
            ma(i)=ma(i)-ma(1+mod(i+30,55))
            if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG
12        continue
13      continue
        inext=0
        inextp=31
        idum=1
      endif
      inext=inext+1
      if(inext.eq.56)inext=1
      inextp=inextp+1
      if(inextp.eq.56)inextp=1
      mj=ma(inext)-ma(inextp)
      if(mj.lt.MZ)mj=mj+MBIG
      ma(inext)=mj
      ran3=mj*FAC
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software m!1.
