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,*)'' write(ifile2,*) '' write(ifile2,*) '
' write(ifile2,*)' ' write(ifile2,*) '' write(ifile2,2179) starname,planetname write(ifile2,*) '' write(ifile2,*) 'quantity | ' write(ifile2,*) 'unit | ' write(ifile2,*) 'value | ' write(ifile2,*) '1-sigma uncertainty | ' write(ifile2,*) '|||
Period | ' write(ifile2,*) 'days | ' write(ifile2,6000) '',Period,' | ' write(ifile2,6000) '',dP,' | ' 6000 format(a16,f12.7,a5) write(ifile2,*) '|||
Time of periastron | ' write(ifile2,*) 'HJD | ' write(ifile2,6001) '',tperi,' | ' write(ifile2,6001) '',dtperi,' | ' 6001 format(a16,f15.7,a5) write(ifile2,*) '|||
Orbital eccentricity | ' write(ifile2,*) 'none | ' write(ifile2,6002) '',e,' | ' write(ifile2,6002) '',de,' | ' 6002 format(a16,f5.3,a5) write(ifile2,*) '|||
Argument of perihelion | ' write(ifile2,*) 'deg | ' write(ifile2,6003) '',360.*om/(2.*pi),' | ' write(ifile2,6003) '',dom,' | ' 6003 format(a16,f7.3,a5) write(ifile2,*) '|||
Velocity half-amplitude | ' write(ifile2,*) 'm s^-1 | ' write(ifile2,6004) '',rk/100.,' | ' write(ifile2,6004) '',drk,' | ' 6004 format(a16,f8.3,a5) write(ifile2,*) '|||
Planetary mass | ' else write(ifile2,*) 'M*sin(i) | ' end if if(rmp.gt.0.3)then write(ifile2,*) 'M_Jup | ' else write(ifile2,*) 'M_Earth | ' endif if(rmp.gt.0.3)then write(ifile2,6005) '',rmp,' | ' else write(ifile2,6005) '',rmp*319.,' | ' endif 6005 format(a16,f6.3,a5) write(ifile2,*) '--- | ' write(ifile2,*) '
Stellar mass | ' write(ifile2,*) 'M_Sun | ' write(ifile2,6006) '',rmstar,' | ' 6006 format(a16,f5.2,a5) write(ifile2,*) '--- | ' write(ifile2,*) '|||
Stellar T_eff | ' write(ifile2,*) 'Kelvin | ' write(ifile2,*) '',int(Teff),' | ' write(ifile2,*) '--- | ' write(ifile2,*) '|||
Planetary T_eff | ' write(ifile2,*) 'Kelvin | ' write(ifile2,*) '',int(Tplan),' | ' write(ifile2,*) '--- | ' write(ifile2,*) '|||
Stellar Radius | ' write(ifile2,*) 'R_Sun | ' write(ifile2,6005) '',rstar,' | ' write(ifile2,*) '--- | ' write(ifile2,*) '|||
Predicted Planetary Radius (w/ +o core) | ' else write(ifile2,*) 'Measured Planetary Radius | +' end if write(ifile2,*) 'R_Jupiter | ' if(itryn.eq.0) then write(ifile2,6005) '',rad(2),' | ' else write(ifile2,6005) '',radknown,' | ' end if write(ifile2,*) '--- | ' write(ifile2,*) '|
T_circ (w/o core) | ' write(ifile2,*) 'Gyr | ' write(ifile2,6005) '',tcirc(2),' | ' write(ifile2,*) '--- | ' write(ifile2,*) '|||
T_damp (w/o core) | ' write(ifile2,*) 'Gyr | ' write(ifile2,6005) '',tdamp(2),' | ' write(ifile2,*) '--- | ' write(ifile2,*) '