program globalupdate
       implicit real*8(a-h,o-z), integer*4(i-n)

c....................................................................

c..    This code provides a master update to the transitsearch.org
c..    transit search project.

c..    This code is undergoing modification as of June 2007

c....................................................................


c.     Physical constants
       parameter(pi=3.14159265358979323846)
       parameter(rmsun=1.98911e+33)
       parameter(rau=1.495978707e+13)
       parameter(G=6.672e-8)
       parameter(sig=5.6703e-5)
       parameter(rlsun=3.826e+33)
       parameter(rmjup=1.8986e+30)
       parameter(radjup=7.1398e+9)
       parameter(rsun=6.696e+10)
       parameter(year=365.25*24.*3600.)

c..    control parmaeters
       parameter(ihbehind=8)

       character*9 starname
       character*9 starkeep
       character*3 cmonth
       character*2 cday
       character*2 cmin
       character*2 chour
       character*1 planetname
       character*1 planetkeep

       character*56 hst1
c      character*21 string1
       character*10 string2
       character*13 string3
       character*13 string9
       character*11 string4
       character*23 string5
       character*23 string10
       character*21 string6
       character*23 string7
       character*21 string8
       character*23 string11
       character*10 string12
       character*3 sstatus
       character*8 string20
       character*18 string21
       character*14 string22 
       character*32 string23
       character*34 string24
       character*57 string25


       character tline1a *4
       character tline2a *15
       character tline2b *5
       character tline3a *24
       character tline3b *2
       character tline3c *9
       character tline4a *15
       character tline4b *5
       character tline5a *15
       character tline5b *5
       character tline6a *15
       character tline6b *5
       character tline7a *15
       character tline7b *5
       character tline8a *15
       character tline8b *5
       character tline9a *16
       character tline9b *5
       character tline10a *15
       character tline10b *5
       character tline11a *24
       character tline11b *20
       character tline12a *24
       character tline12b *18
       character tline13a *5

       dimension rad(4),tcirc(4),tdamp(4)


c..    best guess extrasolar planet properties:

c..    planetary albedos
       Albedo=0.2
c..    tidal quality factors
       Qp=1.0e+6


c..    Julian date for analysis start
       today=real(time()+ihbehind*3600)/86400.
       today=2440587.+today

c..    time (in days) for which prediction tables are made
       wait=2.*(365.25)

c..    open input data files

       open(unit=1,file='orbit.data',status='old')   
       open(unit=2,file='star.data',status='old')
       open(unit=3,file='transit.data',status='old')

c..    open temporary file for output table construction.
       open(unit=12,file='index.blok',status='unknown')



c..    skip header line in orbit.data
       read(1,*) 

c..    read number of active stars
       read(1,*) nstars

c..    skip header line in star.data
       read(2,*)

c..    skip header line in transit.data
       read(3,*)

c..    Begin the main star loop
       do istar=1,nstars

         ifirst=0
         izero=0
         sstatus='out'

c........Read in data from orbit, star, transit dot data files

c..      1. read information from orbit.data file

         read(1,2101) starname,planetname,period,dp,tperi,
     +              dtperi,e,de,om,dom,rk,drk,numplan,itryn

2101     format(a9,a1,1x,f12.7,1x,f12.7,1x,f15.7,1x,f12.7,
     +        1x,f5.3,1x,f5.3,1x,f7.3,1x,f7.3,1x,f6.1,1x,f5.1,1x,i1,1x
     +        ,i1)  


         fsini=1.0
         if(dp.lt.1.0e-6) dp=period*(0.05)


c..      2. read information from the star.data file

         read(2,2105) starname,planetname,rmstar,teff,rstar,
     +              ira1,ira2,idec1,idec2

2105     format(a9,a1,1x,f4.2,1x,f5.0,1x,f5.2,1x,i2,1x,i2,1x,i3,1x,i2)   

c..      3. read information from transit.data for known transits.

         if(itryn.eq.1) then
           read(3,2107) starname,planetname,obsvdur,radknown
         else
           read(3,*) 
         end if
2107     format(a9,a1,1x,f12.7,1x,f7.3)

c....................

c..      set up transit, orbit and results html pages

c         string1='/usr/local/www/laugh/'

         string2=starname(1:9)//planetname(1:1)

         string3='.transits.txt'
         string4='.orbit.html'
         string9='.results.html'
        
         string5=string2(1:10)//string3(1:13)
         string6=string2(1:10)//string4(1:11)
         string10=string2(1:10)//string9(1:13)

         string7=string2(1:10)//string3(1:13)
         string8=string2(1:10)//string4(1:11)
         string11=string2(1:10)//string9(1:13)
       
         string24='cat rblock1 index.blok2 rblock2 > '
         string25=string24(1:34)//string10(1:23)

         ifile=11
         ifile2=61
         ifile3=71

         open(unit=ifile,file=string5,status='unknown')
         open(unit=ifile2,file=string6,status='unknown')
         open(unit=ifile3,file='index.blok2',status='unknown')

c........Completed opening transit, orbit, results pages for current star

c..      find mass, radius estimates, circularization time estimates
c..      semi-major axis:
 
         call circularize(albedo,e,qp,rk,period,rmstar,teff,
     +                    rstar,fsini,rmp,rad,tcirc,tdamp,a,tplan)

c..      compute orbital ephemeris from 
         om=2.*pi*(om/360.)
         f=(pi/2.)-om
         fac=sqrt((1.-e)/(1.+e))
         Ea=2.*atan(fac*tan(f/2.))
         if(Ea.lt.0.) Ea=Ea+2.*pi
         Te=Period*(Ea-e*sin(Ea))/(2.*pi)

c..      now compute transit duration:
         rplansol=rad(2)*(Radjup/rsun)
         df=((Rstar+rplansol)/a)*(rsun/rau)
         ef_fac=sqrt((1-e)/(1+e))
         ef_fac=ef_fac*tan(f/2.)
         ef_fac=2.*atan(ef_fac)
         ef_fac=1-e*cos(ef_fac)
         df=df/ef_fac
         df=asin(df)

         fnew=f+df
         Eanew=2.*atan(fac*tan(fnew/2.))
         if(Eanew.lt.0.) Eanew=Eanew+2.*pi
         Tduration=Period*(Eanew-e*sin(Eanew))/(2.*pi)
         if(tduration.lt.te) tduration=tduration+period
         Tduration=2.*(Tduration-Te)

         if(itryn.eq.1) then
           Tduration=obsvdur/(24.*60.)
         end if


         write(ifile,2175) starname,planetname
2175     format(1x,'Predicted Transit Epochs: ',a9,1x,a1)
         write(ifile,*)   
         if(itryn.eq.0) then

         if ((tduration*24.*60.) .lt. ( 500.)) then
         write(ifile,6009) ' Predicted Duration of a Central Transit: '
     +                    ,tduration*24.*60.,' Minutes'
         else
         write(ifile,6009) ' Predicted Duration of a Central Transit: '
     +                    ,tduration*24.,' Hours'
         end if
6009     format(a42,f6.0,a8)


         else

         if ((tduration*24.*60.) .lt. ( 500.)) then
         write(ifile,6010) 'Observed Transit Duration: '
     +                    ,tduration*24.*60.,' Minutes'
         else
         write(ifile,6010) 'Observed Transit Duration: '
     +                    ,tduration*24.,' Hours'
         end if
6010     format(a27,f6.1,a8)

         end if

         write(ifile,*)   
         write(ifile,3000) 
3000     format(1x,'Begin Transit Window',11x,
     +        'PREDICTED CENTRAL TRANSIT'
     +          ,5x,'End Transit Window')
         write(ifile,3001)
3001     format(37x,'All Times UT')

         write(ifile,*)   
         write(ifile,2176)   
2176     format(31x,'HJD ',7x,'Year',1x,'M',2x,'D',2x,'H',2x,'M')

         do iloop=1,100000
            
c..        central transit
           tjd=Tperi+real(iloop)*period+Te
           iajd=int(tjd+0.5)
           ibjd=int((real(iajd)-1867216.25)/36524.25)
           icjd=iajd+ibjd-int(real(ibjd)/4.)+1525
           idjd=int((real(icjd)-122.1)/365.25)
           iejd=int(365.25*real(idjd))
           ifjd=int((real(icjd)-real(iejd))/30.6001)
           iday=icjd-iejd
           iweird=-int(30.6001*real(ifjd))
           iday=iday+iweird
           day=real(iday)

           fracday=((tjd+0.5)-real(aint(tjd+0.5)))
           hour=real(aint(24.*fracday))
           rminute=60.*((24.*fracday)-real(aint(24.*fracday)))
           rmonth=ifjd-1-12*int(real(ifjd)/14.)
           ryear=idjd-4715-int((7+real(rmonth))/10.)

           if(rmonth.eq.2.and.iday.eq.31) then
             ryear=ryear-1
             iday=iday-2
             day=day-2
           end if
       

c..        early transit

           Tee=(Period-dp)*(Ea-e*sin(Ea))/(2.*pi)-tduration/2.
           tjde=(Tperi-dtperi)+real(iloop)*(period-dp)+Tee
           iajd=int(tjde+0.5)
           ibjd=int((real(iajd)-1867216.25)/36524.25)
           icjd=iajd+ibjd-int(real(ibjd)/4.)+1525
           idjd=int((real(icjd)-122.1)/365.25)
           iejd=int(365.25*real(idjd))
           ifjd=int((real(icjd)-real(iejd))/30.6001)
           iday=icjd-iejd
           iweird=-int(30.6001*real(ifjd))
           iday=iday+iweird
           daye=real(iday)

           fracday=((tjde+0.5)-real(aint(tjde+0.5)))
           houre=real(aint(24.*fracday))
           rminutee=60.*((24.*fracday)-real(aint(24.*fracday)))
           rmonthe=ifjd-1-12*int(real(ifjd)/14.)
           ryeare=idjd-4715-int((7+real(rmonthe))/10.)

           if(rmonthe.eq.2.and.iday.eq.31) then
             ryeare=ryeare-1
             iday=iday-2
             daye=daye-2
           end if

c..        late transit

           Tel=(Period+dp)*(Ea-e*sin(Ea))/(2.*pi)+tduration/2.
           tjdl=(Tperi+dtperi)+real(iloop)*(period+dp)+Tel
           iajd=int(tjdl+0.5)
           ibjd=int((real(iajd)-1867216.25)/36524.25)
           icjd=iajd+ibjd-int(real(ibjd)/4.)+1525
           idjd=int((real(icjd)-122.1)/365.25)
           iejd=int(365.25*real(idjd))
           ifjd=int((real(icjd)-real(iejd))/30.6001)
           iday=icjd-iejd
           iweird=-int(30.6001*real(ifjd))
           iday=iday+iweird
           dayl=real(iday)

           fracday=((tjdl+0.5)-real(aint(tjdl+0.5)))
           hourl=real(aint(24.*fracday))
           rminutel=60.*((24.*fracday)-real(aint(24.*fracday)))
           rmonthl=ifjd-1-12*int(real(ifjd)/14.)
           ryearl=idjd-4715-int((7+real(rmonthl))/10.)

           if(rmonthl.eq.2.and.iday.eq.31) then
             ryearl=ryearl-1
             iday=iday-2
             dayl=dayl-2
           end if

c..        check if we're in transit window:
           if((tjde.lt.today).and.(tjdl.gt.today)) then
            sstatus='in'
           end if

           if( (tjd.gt.(today-120.)) .and. ( ((tjd.lt.today+1000.) 
     +      .and.(period.lt.20.)) .or. ((tjd.lt.today+10000.)               
     +      .and.(period.ge.20.)) ) )then
c..          
             if((ifirst.eq.0).and.(tjd.gt.today)) then
               ifirst=1

c..      month character conversion
         if(int(rmonth).eq.1) cmonth='Jan'
         if(int(rmonth).eq.2) cmonth='Feb'
         if(int(rmonth).eq.3) cmonth='Mar'
         if(int(rmonth).eq.4) cmonth='Apr'
         if(int(rmonth).eq.5) cmonth='May'
         if(int(rmonth).eq.6) cmonth='Jun'
         if(int(rmonth).eq.7) cmonth='Jul'
         if(int(rmonth).eq.8) cmonth='Aug'
         if(int(rmonth).eq.9) cmonth='Sep'
         if(int(rmonth).eq.10) cmonth='Oct'
         if(int(rmonth).eq.11) cmonth='Nov'
         if(int(rmonth).eq.12) cmonth='Dec'

         if(int(day).eq.1) cday='01'
         if(int(day).eq.2) cday='02'
         if(int(day).eq.3) cday='03'
         if(int(day).eq.4) cday='04'
         if(int(day).eq.5) cday='05'
         if(int(day).eq.6) cday='06'
         if(int(day).eq.7) cday='07'
         if(int(day).eq.8) cday='08'
         if(int(day).eq.9) cday='09'
         if(int(day).eq.10) cday='10'
         if(int(day).eq.11) cday='11'
         if(int(day).eq.12) cday='12'
         if(int(day).eq.13) cday='13'
         if(int(day).eq.14) cday='14'
         if(int(day).eq.15) cday='15'
         if(int(day).eq.16) cday='16'
         if(int(day).eq.17) cday='17'
         if(int(day).eq.18) cday='18'
         if(int(day).eq.19) cday='19'
         if(int(day).eq.20) cday='20'
         if(int(day).eq.21) cday='21'
         if(int(day).eq.22) cday='22'
         if(int(day).eq.23) cday='23'
         if(int(day).eq.24) cday='24'
         if(int(day).eq.25) cday='25'
         if(int(day).eq.26) cday='26'
         if(int(day).eq.27) cday='27'
         if(int(day).eq.28) cday='28'
         if(int(day).eq.29) cday='29'
         if(int(day).eq.30) cday='30'
         if(int(day).eq.31) cday='31'

         if(int(hour).eq.0) chour='00'
         if(int(hour).eq.1) chour='01'
         if(int(hour).eq.2) chour='02'
         if(int(hour).eq.3) chour='03'
         if(int(hour).eq.4) chour='04'
         if(int(hour).eq.5) chour='05'
         if(int(hour).eq.6) chour='06'
         if(int(hour).eq.7) chour='07'
         if(int(hour).eq.8) chour='08'
         if(int(hour).eq.9) chour='09'
         if(int(hour).eq.10) chour='10'
         if(int(hour).eq.11) chour='11'
         if(int(hour).eq.12) chour='12'
         if(int(hour).eq.13) chour='13'
         if(int(hour).eq.14) chour='14'
         if(int(hour).eq.15) chour='15'
         if(int(hour).eq.16) chour='16'
         if(int(hour).eq.17) chour='17'
         if(int(hour).eq.18) chour='18'
         if(int(hour).eq.19) chour='19'
         if(int(hour).eq.20) chour='20'
         if(int(hour).eq.21) chour='21'
         if(int(hour).eq.22) chour='22'
         if(int(hour).eq.23) chour='23'
         if(int(hour).eq.24) chour='24'

         if(int(rminute).eq.1) cmin='01'
         if(int(rminute).eq.2) cmin='02'
         if(int(rminute).eq.3) cmin='03'
         if(int(rminute).eq.4) cmin='04'
         if(int(rminute).eq.5) cmin='05'
         if(int(rminute).eq.6) cmin='06'
         if(int(rminute).eq.7) cmin='07'
         if(int(rminute).eq.8) cmin='08'
         if(int(rminute).eq.9) cmin='09'
         if(int(rminute).eq.10) cmin='11'
         if(int(rminute).eq.11) cmin='11'
         if(int(rminute).eq.12) cmin='12'
         if(int(rminute).eq.13) cmin='13'
         if(int(rminute).eq.14) cmin='14'
         if(int(rminute).eq.15) cmin='15'
         if(int(rminute).eq.16) cmin='16'
         if(int(rminute).eq.17) cmin='17'
         if(int(rminute).eq.18) cmin='18'
         if(int(rminute).eq.19) cmin='19'
         if(int(rminute).eq.20) cmin='20'
         if(int(rminute).eq.21) cmin='21'
         if(int(rminute).eq.22) cmin='22'
         if(int(rminute).eq.23) cmin='23'
         if(int(rminute).eq.24) cmin='24'
         if(int(rminute).eq.25) cmin='25'
         if(int(rminute).eq.26) cmin='26'
         if(int(rminute).eq.27) cmin='27'
         if(int(rminute).eq.28) cmin='28'
         if(int(rminute).eq.29) cmin='29'
         if(int(rminute).eq.30) cmin='30'
         if(int(rminute).eq.31) cmin='31'
         if(int(rminute).eq.32) cmin='32'
         if(int(rminute).eq.33) cmin='33'
         if(int(rminute).eq.34) cmin='34'
         if(int(rminute).eq.35) cmin='35'
         if(int(rminute).eq.36) cmin='36'
         if(int(rminute).eq.37) cmin='37'
         if(int(rminute).eq.38) cmin='38'
         if(int(rminute).eq.39) cmin='39'
         if(int(rminute).eq.40) cmin='40'
         if(int(rminute).eq.41) cmin='41'
         if(int(rminute).eq.42) cmin='42'
         if(int(rminute).eq.43) cmin='43'
         if(int(rminute).eq.44) cmin='44'
         if(int(rminute).eq.45) cmin='45'
         if(int(rminute).eq.46) cmin='46'
         if(int(rminute).eq.47) cmin='47'
         if(int(rminute).eq.48) cmin='48'
         if(int(rminute).eq.49) cmin='49'
         if(int(rminute).eq.50) cmin='50'
         if(int(rminute).eq.51) cmin='51'
         if(int(rminute).eq.52) cmin='52'
         if(int(rminute).eq.53) cmin='53'
         if(int(rminute).eq.54) cmin='54'
         if(int(rminute).eq.55) cmin='55'
         if(int(rminute).eq.56) cmin='56'
         if(int(rminute).eq.57) cmin='57'
         if(int(rminute).eq.58) cmin='58'
         if(int(rminute).eq.59) cmin='59'
         if(int(rminute).eq.60) cmin='60'


2103       format(a2,':',a2,1x,a3,1x,a2,',',1x,i4)
           isaveyear=ryear
 
             endif 

             write(ifile,2177) tjde,
     +       int(ryeare),int(rmonthe),int(daye),int(houre),int(rminutee)
     +      ,tjd,int(ryear),int(rmonth),int(day),
     +       int(hour),int(rminute)
     +      ,tjdl,int(ryearl),int(rmonthl),int(dayl),
     +       int(hourl),int(rminutel)

2177         format(1x,f10.2,1x,i4,1x,i2,1x,i2,1x,i2,1x,i2,3x,
     +                 f10.2,1x,i4,1x,i2,1x,i2,1x,i2,1x,i2,3x,
     +                 f10.2,1x,i4,1x,i2,1x,i2,1x,i2,1x,i2)

           end if

            if((tjd.gt.today+730.  ).and.(period.lt.20.)) goto 2102
            if((tjd.gt.today+7300. ).and.(period.gt.20.)) goto 2102

         end do

         

2102     continue



c..      planetary data output block
         write(ifile2,*)'<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 T
     +ransitional//EN">'
         write(ifile2,*) '<html>'
         write(ifile2,*) '<head>'
         write(ifile2,*)' <meta http-equiv="content-type" content="text
     +/html;charset=ISO-8859-1">'
         write(ifile2,*) '<meta name="generator" content="predict.f">'
         write(ifile2,2179) starname,planetname
         write(ifile2,*) '</head>'

         write(ifile2,*) '<br>'
         write(ifile2,2178) starname,planetname
         write(ifile2,*) '<br>'
         write(ifile2,*) '<br>'
2178     format('Orbital Data and Parent Star Data for: ',a9,1x,a1)
2179     format('<title>Data for ',a9,1x,a1,'</title>')

         write(ifile2,*) '<body bgcolor="#ffffff">'
         write(ifile2,*) '<table width="637" border="1" cellspacing="2" 
     +cellpadding="0">'
         write(ifile2,*) '<tr>'
         write(ifile2,*) '<td width="261"><b>quantity</b></td>'
         write(ifile2,*) '<td width="69"><b>unit</b></td>'
         write(ifile2,*) '<td width="135"><b>value</b></td>'
         write(ifile2,*) '<td width="150"><b>1-sigma uncertainty</b></td
     +>'
         write(ifile2,*) '</tr>'

         write(ifile2,*) '<tr>'
         write(ifile2,*) '<td width="261">Period</td>'
         write(ifile2,*) '<td width="69">days</td>'
         write(ifile2,6000) '<td width="135">',Period,'</td>'
         write(ifile2,6000) '<td width="150">',dP,'</td>'
6000     format(a16,f12.7,a5)
         write(ifile2,*) '</tr>'
         write(ifile2,*) '<tr>'
         write(ifile2,*) '<td width="261">Time of periastron</td>'
         write(ifile2,*) '<td width="69">HJD</td>'
         write(ifile2,6001) '<td width="135">',tperi,'</td>'
         write(ifile2,6001) '<td width="150">',dtperi,'</td>'
6001     format(a16,f15.7,a5)
         write(ifile2,*) '</tr>'
         write(ifile2,*) '<tr>'
         write(ifile2,*) '<td width="261">Orbital eccentricity</td>'
         write(ifile2,*) '<td width="69">none</td>'
         write(ifile2,6002) '<td width="135">',e,'</td>'
         write(ifile2,6002) '<td width="150">',de,'</td>'
6002     format(a16,f5.3,a5)
         write(ifile2,*) '</tr>'
         write(ifile2,*) '<tr>'
         write(ifile2,*) '<td width="261">Argument of perihelion</td>'
         write(ifile2,*) '<td width="69">deg</td>'
         write(ifile2,6003) '<td width="135">',360.*om/(2.*pi),'</td>'
         write(ifile2,6003) '<td width="150">',dom,'</td>'
6003     format(a16,f7.3,a5)
         write(ifile2,*) '</tr>'
         write(ifile2,*) '<tr>'
         write(ifile2,*) '<td width="261">Velocity half-amplitude</td>'
         write(ifile2,*) '<td width="69">m s^-1</td>'
         write(ifile2,6004) '<td width="135">',rk/100.,'</td>'
         write(ifile2,6004) '<td width="150">',drk,'</td>'
6004     format(a16,f8.3,a5)
         write(ifile2,*) '</tr>'
         write(ifile2,*) '<tr>'

         if(itryn.eq.1) then
           write(ifile2,*) '<td width="261">Planetary mass</td>'
         else
           write(ifile2,*) '<td width="261">M*sin(i)</td>'
         end if

         if(rmp.gt.0.3)then
           write(ifile2,*) '<td width="69">M_Jup</td>'
         else
           write(ifile2,*) '<td width="69">M_Earth</td>'
         endif

         if(rmp.gt.0.3)then
           write(ifile2,6005) '<td width="135">',rmp,'</td>'
         else
           write(ifile2,6005) '<td width="135">',rmp*319.,'</td>'
         endif
6005     format(a16,f6.3,a5)

         write(ifile2,*) '<td width="150">---</td>'
         write(ifile2,*) '</tr>'
         write(ifile2,*) '<tr>'
         write(ifile2,*) '<td width="261">Stellar mass</td>'
         write(ifile2,*) '<td width="69">M_Sun</td>'
         write(ifile2,6006) '<td width="135">',rmstar,'</td>'
6006     format(a16,f5.2,a5)
         write(ifile2,*) '<td width="150">---</td>'
         write(ifile2,*) '</tr>'
         write(ifile2,*) '<tr>'
         write(ifile2,*) '<td width="261">Stellar T_eff</td>'
         write(ifile2,*) '<td width="69">Kelvin</td>'
         write(ifile2,*) '<td width="135">',int(Teff),'</td>'
         write(ifile2,*) '<td width="150">---</td>'
         write(ifile2,*) '</tr>'
         write(ifile2,*) '<tr>'
         write(ifile2,*) '<td width="261">Planetary T_eff</td>'
         write(ifile2,*) '<td width="69">Kelvin</td>'
         write(ifile2,*) '<td width="135">',int(Tplan),'</td>'
         write(ifile2,*) '<td width="150">---</td>'
         write(ifile2,*) '</tr>'
         write(ifile2,*) '<tr>'
         write(ifile2,*) '<td width="261">Stellar Radius</td>'
         write(ifile2,*) '<td width="69">R_Sun</td>'
         write(ifile2,6005) '<td width="135">',rstar,'</td>'
         write(ifile2,*) '<td width="150">---</td>'
         write(ifile2,*) '</tr>'
         write(ifile2,*) '<tr>'
         if(itryn.eq.0) then
         write(ifile2,*) '<td width="261">Predicted Planetary Radius (w/
     +o core)</td>'
         else
         write(ifile2,*) '<td width="261">Measured Planetary Radius</td>
     +'
         end if
         write(ifile2,*) '<td width="69">R_Jupiter</td>'
         if(itryn.eq.0) then
           write(ifile2,6005) '<td width="135">',rad(2),'</td>'
         else
           write(ifile2,6005) '<td width="135">',radknown,'</td>'
         end if
         write(ifile2,*) '<td width="150">---</td>'
         write(ifile2,*) '</tr>'
         write(ifile2,*) '<tr>'
         write(ifile2,*) '<td width="261">T_circ (w/o core)</td>'
         write(ifile2,*) '<td width="69">Gyr</td>'
         write(ifile2,6005) '<td width="135">',tcirc(2),'</td>'
         write(ifile2,*) '<td width="150">---</td>'
         write(ifile2,*) '</tr>'
         write(ifile2,*) '<tr>'
         write(ifile2,*) '<td width="261">T_damp (w/o core)</td>'
         write(ifile2,*) '<td width="69">Gyr</td>'
         write(ifile2,6005) '<td width="135">',tdamp(2),'</td>'
         write(ifile2,*) '<td width="150">---</td>'
         write(ifile2,*) '</tr>'

         write(ifile2,*) '</table>'



         write(ifile2,*) '<br>'
         write(ifile2,*) '<br>'

         write(ifile2,*) 'Planet radii from theoretical models of:'
         write(ifile2,*) '<br>'
         write(ifile2,*) 'Bodenheimer, Laughlin, and Lin 2003, Astrophys
     +ical Journal, v. 592, 555-563.'
         write(ifile2,*) '<br>'
         write(ifile2,*) '<br>'

         write(ifile2,*) 'Damping timescales from:'
         write(ifile2,*) '<br>'
         hst1='<a href="http://arxiv.org/abs/astro-ph/0606346">link</a>'   
         write(ifile2,*) 'Adams and Laughlin 2006, ApJ : ',hst1
         write(ifile2,*) '<br>'
         write(ifile2,*) '<br>'
         write(ifile2,*) '</p>'
         write(ifile2,*) '</body>'
         write(ifile2,*) '<script'
         write(ifile2,*) 'src='
         write(ifile2,*) '"http://www.google-analytics.com/urchin.js"'  
         write(ifile2,*) 'type="text/javascript">'
         write(ifile2,*) '</script>'
         write(ifile2,*) '<script type="text/javascript">'
         write(ifile2,*) '_uacct = "UA-2446040-4";'
         write(ifile2,*) 'urchinTracker();'
         write(ifile2,*) '</script>'
         write(ifile2,*) '</body>'
         write(ifile2,*) '</html>'

c....    end data webpage

c..      compute transit probability
         Prob=0.0045*(1./a)*((Rstar*Rsun+Rad(2)*Radjup)/Rsun)
         Prob=Prob*(1+e*cos((pi/2.)-om))/(1-e**2)

c..      check if transiting
         if(itryn.eq.1) then
           prob=100.
         else
           probtot=probtot+prob
           prob=100.*prob
         end if
         

c..      compute central transit depth
         depth=((Rad(2)*Radjup)**2)/((Rstar*Rsun)**2)
         depth=depth*100.
         if(itryn.eq.1) then
           depth=((Radknown*Radjup)**2)/((Rstar*Rsun)**2)
           depth=depth*100.
           Rad(2)=Radknown
         end if

c..      write html code block here:
         tline1a='<tr>'
         write(12,4001) tline1a
4001     format(a4)

         tline2a='<td width="91">'
         tline2b='</td>'
         write(12,4002) tline2a,starname,tline2b
4002     format(a15,a9,a5)

         tline3a='<td width="50"><a href="'
         tline3b='">'
         tline3c='</a></td>'
         write(12,4003) tline3a,string8,tline3b,planetname,tline3c
4003     format(a24,a21,a2,a1,a9)

         tline4a='<td width="6%">'
         tline4b='</td>'
         if(period.lt.10.) then
           write(12,4004) tline4a,period,tline4b
4004       format(a15,f5.3,a5)
         elseif(period.ge.10. .and. period.lt.100.) then
           write(12,4005) tline4a,period,tline4b
4005       format(a15,f5.2,a5)
         elseif(period.ge.100. .and. period.lt.1000.) then
           write(12,4006) tline4a,period,tline4b
4006       format(a15,f5.1,a5)
         elseif(period.ge.1000.) then
           write(12,4007) tline4a,period,tline4b
4007       format(a15,f5.0,a5)
         end if

         tline5a='<td width="6%">'
         tline5b='</td>'
         write(12,4008) tline5a,prob,tline5b
4008     format(a15,f5.1,a5)

         tline6a='<td width="55">'
         tline6b='</td>'

         if(ira1.ge.10.and.ira2.ge.10) then
           write(12,4009) tline6a,ira1,ira2,tline6b
4009       format(a15,i2,':',i2,a5)
         elseif(ira1.ge.10.and.ira2.lt.10) then
           write(12,4109) tline6a,ira1,ira2,tline6b
4109       format(a15,i2,':0',i1,a5)
         elseif(ira1.lt.10.and.ira2.lt.10) then
           write(12,4209) tline6a,ira1,ira2,tline6b
4209       format(a15,'0',i1,':0',i1,a5)
         else
           write(12,4309) tline6a,ira1,ira2,tline6b
4309       format(a15,'0',i1,':',i2,a5)
         end if

        
         

         tline7a='<td width="55">'
         tline7b='</td>'

         if(idec1.ge.0.) then

           if(idec1.ge.10.and.idec2.ge.10) then
             write(12,4010) tline7a,idec1,idec2,tline7b
4010         format(a15,'+',i2,':',i2,a5)
           elseif(idec1.ge.10.and.idec2.lt.10) then
             write(12,4110) tline7a,idec1,idec2,tline7b
4110         format(a15,'+',i2,':0',i1,a5)
           elseif(idec1.lt.10.and.idec2.lt.10) then
             write(12,4210) tline7a,idec1,idec2,tline7b
4210         format(a15,'+0',i1,':0',i1,a5)
           else
             write(12,4310) tline7a,idec1,idec2,tline7b
4310         format(a15,'+0',i1,':',i2,a5)
           end if
    
         else

           idec1=abs(idec1)
           if(idec1.ge.10.and.idec2.ge.10) then
             write(12,4410) tline7a,idec1,idec2,tline7b
4410         format(a15,'-',i2,':',i2,a5)
           elseif(idec1.ge.10.and.idec2.lt.10) then
             write(12,4510) tline7a,idec1,idec2,tline7b
4510         format(a15,'-',i2,':0',i1,a5)
           elseif(idec1.lt.10.and.idec2.lt.10) then
             write(12,4610) tline7a,idec1,idec2,tline7b
4610         format(a15,'-0',i1,':0',i1,a5)
           else
             write(12,4710) tline7a,idec1,idec2,tline7b
4710         format(a15,'-0',i1,':',i2,a5)
           end if
   
         end if

         tline8a='<td width="38">'
         tline8b='</td>'
         if(depth.ge.0.1) then
           write(12,4011) tline8a,depth,tline8b
4011       format(a15,f5.2,a5)
         else
           write(12,4111) tline8a,depth,tline8b
4111       format(a15,f5.3,a5)
         end if

         tline9a='<td width="202">'
         tline9b='</td>'
         write(12,4012) tline9a,chour,cmin,cmonth,cday,isaveyear,tline9b          
4012     format(a16,a2,':',a2,1x,a3,1x,a2,',',1x,i4,a5)

         tline10a='<td width="60">'
         tline10b='</td>'
         write(12,4013) tline10a,sstatus,tline10b
4013     format(a15,a3,a5)

         tline11a='<td width="83"><a href="'
         tline11b='">Ephemeris</a></td>'
         write(12,4014) tline11a,string7,tline11b
4014     format(a24,a23,a20)

         tline12a='<td width="91"><a href="'
         tline12b='">Results</a></td>'
         write(12,4015) tline12a,string11,tline12b
4015     format(a24,a23,a18)

         tline13a='</tr>'
         write(12,4016) tline13a
4016     format(a5)

       close(ifile)
       close(ifile2)

c..    Write the default results html file
      if(itryn.eq.0) then
      write(ifile3,*)
     +'A results page has not yet been created for this star.'    
      else
      write(ifile3,*)
     +'A results page has not yet been created for this known transiting      
     + planet.'
      endif
      close(ifile3)
       call system(string25)  

c..    end of star loop
       end do

c..    close up the table after the loop finished
       close(12)


c..    rebuild the candidates table
       call system('cat block1 index.blok block2 > candidates.html')  

c..    reinstate existing results files
       call system('./resultscopy')

       end
     
  


       subroutine circularize(Albedo,ecc,Qp,rk,period,rmstar,tstar,
     +    rstar,fsini,rmp,rad,tcirc,tdamp,a,tplan)

       implicit real*8(a-h,o-z), integer*4(i-n)
       dimension tcirc(4),tdamp(4),rad(4)

       parameter(pi=3.14159265358979323846)

       rmsun=1.98911e+33
       rau=1.495978707e+13
       G=6.672e-8
       sig=5.6703e-5
       rLsun=3.826e+33
       rmjup=1.8986e+30
       radjup=7.1398e+9
       rsun=6.696e+10
       year=365.25*24.*3600.

c..    convert quantities to cgs
       rmstar=rmstar*rmsun
       rstar=rstar*rsun
       period=period*24.*3600.
       rn=2.*pi/(period)
       rK=rK*100.

c..    iterate to get planet m sin(i) from radial velocity K
       do i=1,100000000
         rmp=(real(i)/10000.)
         rmp=rmp*rmjup

         abod=period*period*G*(rmstar+rmp)
         abod=abod/(4.*pi*pi)
         abod=abod**0.33333333333
         q1=(1./rK)*(2.*pi*G/period)**0.333333333
         q1=q1/sqrt(1-ecc*ecc)
         q1=(q1*rmp*fsini)/(rmstar+rmp)**0.66666666
         if(q1.gt.1.) then
           goto 4104
         end if
       end do
4104   continue


  
c..    semi-major axis
       a=(G*(rmp+rmstar)*(period**2)/(4.*pi**2))**0.333333


c..    compute effective temperature, Teff for the planet
       rlstar=4.*pi*rstar**2*sig*(tstar**4)
       Teff=(1.-Albedo)*rLstar
       fac=(a**2)*(1-(ecc**2))**0.5
       Teff=Teff/(sig*16.*pi*fac)
       Teff=Teff**0.25
       tplan=teff

c..    compute radius, circularization timescales
c..    code: iflag=1 core, no kinetic heating
c..    code: iflag=2 no core, no kinetic heating
c..    code: iflag=3 core, kinetic heating
c..    code: iflag=4 no core, kinetic heating

       do iflag=1,4
         call radcompute(teff,rmp/rmjup,radplanet,iflag)
         radplanet=radplanet*radjup

c..      tidal circularization timescales:

c..      first expression
         tep=((4.*Qp)/(63.*rn))*(rmp/rmstar)*(a/radplanet)**5
         tep=tep*(1-ecc*2)**(13./2.)*(1./(1+6.*ecc**2))
         tep=(tep/year)/(10**9)

c..      second expression
c         tep2=0.63*(rmp/rmjup)*(rmsun/rmstar)**1.5*
c     +      (a/(10.*rsun))**6.5*(radjup/radplanet)**5

         rad(iflag)=radplanet/radjup
c         tcirc(iflag)=(tep+tep2)/2.
         tcirc(iflag)=tep
         tdamp(iflag)=((1-ecc**2)/(2*ecc**2))*tep

       end do




       rmp=rmp/rmjup

       if (rmp.lt.0.025) then
         do iflag=1,4
           rad(iflag)=((0.75*rmp*rmjup/(pi*5.))**0.333333)/radjup
         end do
       end if

       a=a/rau 
       period=period/(24.*3600.)
       rmstar=rmstar/rmsun
       rstar=rstar/rsun

       end
     
  




           Subroutine radcompute(T,rho,k,iflag)

c..........This subroutine computes the radius for a specificed temperature
c..........and mass by linear interpolation between the proper points
c..........of the radius table in common.

           implicit real*8(a-h,o-z), integer*4(i-n)
           real*8 k,krare,kdense

           dimension rhoref(6),tref(5),rkapparef(4,6,5)
           data irhs,its /6,5/
           
           data rhoref /0.10,0.23,0.69,1.00,3.00,25.00/
           data tref /85,500,1000,1500,2000/

c..        m=0.11
           rkapparef(1,1,1)=0.61
           rkapparef(1,1,2)=0.66
           rkapparef(1,1,3)=0.69
           rkapparef(1,1,4)=0.74
           rkapparef(1,1,5)=0.87

           rkapparef(2,1,1)=0.89
           rkapparef(2,1,2)=1.01
           rkapparef(2,1,3)=1.09
           rkapparef(2,1,4)=1.20
           rkapparef(2,1,5)=1.75

           rkapparef(3,1,1)=0.61
           rkapparef(3,1,2)=0.80
           rkapparef(3,1,3)=1.14
           rkapparef(3,1,4)=1.69
           rkapparef(3,1,5)=2.50

           rkapparef(4,1,1)=0.89
           rkapparef(4,1,2)=1.13
           rkapparef(4,1,3)=1.72
           rkapparef(4,1,4)=2.50
           rkapparef(4,1,5)=3.50

c..        m=0.23
           rkapparef(1,2,1)=0.82
           rkapparef(1,2,2)=0.88
           rkapparef(1,2,3)=0.90
           rkapparef(1,2,4)=0.95
           rkapparef(1,2,5)=1.07

           rkapparef(2,2,1)=0.95
           rkapparef(2,2,2)=1.03
           rkapparef(2,2,3)=1.07
           rkapparef(2,2,4)=1.14
           rkapparef(2,2,5)=1.31

           rkapparef(3,2,1)=0.82
           rkapparef(3,2,2)=1.02
           rkapparef(3,2,3)=1.34
           rkapparef(3,2,4)=1.61
           rkapparef(3,2,5)=2.50

           rkapparef(4,2,1)=0.95
           rkapparef(4,2,2)=1.14
           rkapparef(4,2,3)=1.53
           rkapparef(4,2,4)=1.80
           rkapparef(4,2,5)=2.50

c..        m=0.69
           rkapparef(1,3,1)=0.95
           rkapparef(1,3,2)=1.00
           rkapparef(1,3,3)=1.01
           rkapparef(1,3,4)=1.03
           rkapparef(1,3,5)=1.10

           rkapparef(2,3,1)=1.03
           rkapparef(2,3,2)=1.09
           rkapparef(2,3,3)=1.10
           rkapparef(2,3,4)=1.13
           rkapparef(2,3,5)=1.22

           rkapparef(3,3,1)=0.95
           rkapparef(3,3,2)=1.02
           rkapparef(3,3,3)=1.16
           rkapparef(3,3,4)=1.35
           rkapparef(3,3,5)=1.74

           rkapparef(4,3,1)=1.03
           rkapparef(4,3,2)=1.11
           rkapparef(4,3,3)=1.27
           rkapparef(4,3,4)=1.51
           rkapparef(4,3,5)=1.81

c..        m=1.00
           rkapparef(1,4,1)=1.01
           rkapparef(1,4,2)=1.05
           rkapparef(1,4,3)=1.06
           rkapparef(1,4,4)=1.08
           rkapparef(1,4,5)=1.14

           rkapparef(2,4,1)=1.05
           rkapparef(2,4,2)=1.10
           rkapparef(2,4,3)=1.11
           rkapparef(2,4,4)=1.13
           rkapparef(2,4,5)=1.20

           rkapparef(3,4,1)=1.01
           rkapparef(3,4,2)=1.06
           rkapparef(3,4,3)=1.13
           rkapparef(3,4,4)=1.30
           rkapparef(3,4,5)=1.47

           rkapparef(4,4,1)=1.05
           rkapparef(4,4,2)=1.11
           rkapparef(4,4,3)=1.18
           rkapparef(4,4,4)=1.38
           rkapparef(4,4,5)=1.61

c..        m=3.00
           rkapparef(1,5,1)=1.07
           rkapparef(1,5,2)=1.10
           rkapparef(1,5,3)=1.11
           rkapparef(1,5,4)=1.12
           rkapparef(1,5,5)=1.18

           rkapparef(2,5,1)=1.07
           rkapparef(2,5,2)=1.11
           rkapparef(2,5,3)=1.11
           rkapparef(2,5,4)=1.13
           rkapparef(2,5,5)=1.19

           rkapparef(3,5,1)=1.07
           rkapparef(3,5,2)=1.10
           rkapparef(3,5,3)=1.13
           rkapparef(3,5,4)=1.25
           rkapparef(3,5,5)=1.40

           rkapparef(4,5,1)=1.07
           rkapparef(4,5,2)=1.11
           rkapparef(4,5,3)=1.14
           rkapparef(4,5,4)=1.26
           rkapparef(4,5,5)=1.35

c..        m=15.00
           rkapparef(1,6,1)=1.07
           rkapparef(1,6,2)=1.10
           rkapparef(1,6,3)=1.11
           rkapparef(1,6,4)=1.12
           rkapparef(1,6,5)=1.18

           rkapparef(2,6,1)=1.07
           rkapparef(2,6,2)=1.11
           rkapparef(2,6,3)=1.11
           rkapparef(2,6,4)=1.13
           rkapparef(2,6,5)=1.19

           rkapparef(3,6,1)=1.07
           rkapparef(3,6,2)=1.10
           rkapparef(3,6,3)=1.13
           rkapparef(3,6,4)=1.25
           rkapparef(3,6,5)=1.40

           rkapparef(4,6,1)=1.07
           rkapparef(4,6,2)=1.11
           rkapparef(4,6,3)=1.14
           rkapparef(4,6,4)=1.26
           rkapparef(4,6,5)=1.35

2109  format(19f7.3)
           if ((T.le.Tref(1)).or.(T.ge.Tref(its)).or.(Rho.le.Rhoref(1))
     +     .or.(Rho.ge.Rhoref(irhs))) then
c..           default value for small planets
c              k=0.35
              k=0.1*(rho*319.)**0.3333
              goto 2104
           endif

c..........Determine which of the reference points bracket the temperature.

           do i=1,its
               if (Tref(i).gt.T) then
                   ihot = i
                   icool = i-1
                   goto 10
               endif
           end do

10         continue

           do j=1,irhs
               if (Rho.lt.Rhoref(j)) then
                   jrare = j-1
                   jdense = j
                   goto 20
               end if
           end do

20         continue

c..........Determine Krare and Kdense by a linear interpolation in temperature.

        krare = rkapparef(iflag,jrare,icool)
     +  +((rkapparef(iflag,jrare,ihot)-rkapparef(iflag,jrare,icool))/
     +  (Tref(ihot)-Tref(icool)))*(T-Tref(icool))

        kdense = rkapparef (iflag,jdense,icool)
     +  +((rkapparef(iflag,jdense,ihot)-rkapparef(iflag,jdense,icool))/
     +  (Tref(ihot)-Tref(icool)))*(T-Tref(icool))

c.........Determine kappa by a linear interpolation in the density

        k = krare + ((kdense-krare)/(rhoref(jdense)-rhoref(jrare)))
     +  *(rho - rhoref(jrare))

c.........We have thus determined k.
 
2104   continue

       return
       end