program solar_eclipse implicit none save c..exercises the solar eclipse program character*15 what(6) character*3 amon(12) integer iphase,iyear,imonth,iday,nstep, 1 iy,im,id,i double precision tjd,lambda,phi,tumbra,rhour,tstep,tjdin, 1 rh,xh,xm,xs,phid,phim,phis,lamd,lamm,lams, 2 tumh,tumm,tums,tjd_beg,tjd_end,tday c..for picking which ephemeris to use integer which_eph common /dop/ which_eph c..for easy phase identification integer noeclipse,partial,non_cen_ann,non_cen_tot, 1 annular,total parameter (noeclipse = 1, 1 partial = 2, 2 non_cen_ann = 3, 3 non_cen_tot = 4, 4 annular = 5, 5 total = 6) c..initialize data what/'noeclipse', 'partial', 'non-cen annular', 1 'non-cen total', 'annular', 'total' / data amon /'jan' , 'feb' , 'mar' , 'apr' , 'may' , 'jun' , 1 'jul' , 'aug' , 'sep' , 'oct' , 'nov' , 'dec' / c..popular format statements 01 format(1x,t2,a,t12,a,t20,a,t30,a,t41,a,t52,a) 02 format(1x,i2.2,a,i4.4,' ',i2.2,i3.2, 1 ' ',i3.2,i4.2,' ',i4.2,i4.2, 2 ' ',0p1f5.2,' ',a) c..initialize the name database which_eph = 406 call init_names(0) c..get the input write(6,*) 'give year month day hour of new moon', 1 ' and tstep in minutes nstep =>' read(5,*) iyear,imonth,iday,rhour,tstep c..convert the calander date of the new moon into a ut julian date call juldat(iyear,imonth,iday,rhour,tjd) c..write a header write(6,*) write(6,01) 'date','time','latitude','longitude','duration', 1 'phase' write(6,01) ' ','h m','degrees','degrees','minutes',' ' c..figure the start end end times tday = tstep/1440.0d0 c tjd_beg = tjd tjd_beg = tjd - 0.13d0 tjd_end = tjd + 0.13d0 nstep = (tjd_end - tjd_beg)/tday + 1 c..for each time point do i=1,nstep tjdin = tjd_beg + float(i-1)*tday c..the the geographical coordinates and duration time of the centerline call central(tjdin,lambda,phi,tumbra,iphase) c..say what we got call caldat(tjdin,iy,im,id,rh) call dh2h(rh,xh,xm,xs) if (xs .gt. 30.0) xm = xm + 1.0d0 call dh2h(phi,phid,phim,phis) if (phis .gt. 30.0) phim = phim + 1.0d0 call dh2h(lambda,lamd,lamm,lams) if (lams .gt. 30.0) lamm = lamm + 1.0d0 write(6,02) id,amon(im),iyear,int(xh),int(xm), 1 int(phid),int(phim),int(lamd),int(lamm), 2 tumbra,what(iphase) enddo end subroutine central(tjd,lambda,phi,tumbra,iphase) implicit none save c..computes the central line, phase and duration of a solar eclipse c..input: c..tjd = ut julian date of interst c..output: c..lambda = geographic longitude of the shadow center in degrees c..phi = geographic latitude of the shadow center in degrees c..tumbra = duration of the total or annular phase at centerline in min c..iphase = phase of the eclipse c.. 1 = no eclipse c.. 2 = partial c.. 3 = non-central annular c.. 4 = non-central total c.. 5 = annular c.. 6 = total c..declare the pass integer iphase double precision tjd,lambda,phi,tumbra c..for easy phase identification integer noeclipse,partial,non_cen_ann,non_cen_tot, 1 annular,total parameter (noeclipse = 1, 1 partial = 2, 2 non_cen_ann = 3, 3 non_cen_tot = 4, 4 annular = 5, 5 total = 6) c..a small time interval in minutes double precision dt parameter (dt = 0.1d0) c..number of minutes in a julian day double precision minjul parameter (minjul = 24.0d0 * 60.0d0) c..angular velocity of earth in rad/min double precision omega parameter (omega = 4.3755d-3) c..local variables integer ip double precision xs,ys,zs,rs,xm,ym,zm,rm, 1 xshad,yshad,zshad,ex,ey,ez,dumbra, 2 tjdh,tjdl,gst,lst,pos(3),tjd1,ra,dec, 3 w,xx,yy,zz,exx,eyy,ezz,dx,dy,dz,dd,du c..popular format statements 222 format(1x,1p6e14.6) c..get the position of the sun tjd1 = tjd call sunmoon(tjd1,xm,ym,zm,rm,xs,ys,zs,rs) c..local mean sidereal time in hours at zero longitude tjdh = int(tjd1) tjdl = tjd1 - tjdh call sidtim(tjdh,tjdl,0,gst) lst = gst c..find the intersection of the shadow axis with the earth call intersect(xm,ym,zm,rm,xs,ys,zs,rs, 1 xshad,yshad,zshad,ex,ey,ez,dumbra,iphase) c..for central phases only if (iphase .lt. annular) then lambda = 0.0d0 phi = 0.0d0 tumbra = 0.0d0 else c..get the geographic coordinates pos(1) = xshad pos(2) = yshad pos(3) = zshad call angles(pos,ra,dec) phi = dec + 0.1924d0 * sind(2.0d0*dec) lambda = 15.0d0 * (lst - ra) if (lambda .gt. 180.0) lambda = lambda - 360.0d0 if (lambda .lt. -180.0) lambda = lambda + 360.0d0 c..get the duration for this location c..get shadow coordinates at t+dt w = dt * omega tjd1 = tjd + dt/minjul call sunmoon(tjd1,xm,ym,zm,rm,xs,ys,zs,rs) call intersect(xm,ym,zm,rm,xs,ys,zs,rs, 1 xx,yy,zz,exx,eyy,ezz,du,ip) c..if the positive delta went overboard, do the negative if (ip .le. annular) then w = -dt * omega tjd1 = tjd - dt/minjul call sunmoon(tjd1,xm,ym,zm,rm,xs,ys,zs,rs) call intersect(xm,ym,zm,rm,xs,ys,zs,rs, 1 xx,yy,zz,exx,eyy,ezz,du,ip) end if c..displacement vector of the shadow on earth dx = xx - xshad + w*yshad dy = yy - yshad - w*xshad dz = zz - zshad dd = sqrt(dx**2 + dy**2 + dz**2 - (dx*ex + dy*ey + dz*ez)**2) tumbra = dt * abs(dumbra)/dd end if return end subroutine intersect(xm1,ym1,zm1,rm1,xs1,ys1,zs1,rs1, 1 xshad,yshad,zshad,ex,ey,ez,dumbra,iphase) implicit none save c..calculates the intersection of the moons shadow axis with the earth c..input: c..xm1 ym1 zm1 = geocentric equatorial position vector of moon in au c..rm1 = geocentic distance to moon in au c..xs1 ys1 zs1 = geocentric equatorial position vector of sun in au c..rs1 = geocentic distance to sun in au c..output: c..xshad yshad zshad = geocentric equatorial position vector of shadow point c.. in equatorial earth radii c..ex ey ez = unit vector of the shadow axis c..dumbra = umbra diameter in equatorial earth radii c..iphase = type of event c.. 1 = no eclipse c.. 2 = partial c.. 3 = non-central annular c.. 4 = non-central total c.. 5 = annular c.. 6 = total c..declare the pass integer iphase double precision xm1,ym1,zm1,rm1,xs1,ys1,zs1,rs1, 1 xshad,yshad,zshad,ex,ey,ez,dumbra c..local variables double precision xm,ym,zm,rm,xs,ys,zs,rs, 1 xms,yms,zms,dms,s0,delta,r0, 2 s,dpenumbra c..for easy phase identification integer noeclipse,partial,non_cen_ann,non_cen_tot, 1 annular,total parameter (noeclipse = 1, 1 partial = 2, 2 non_cen_ann = 3, 3 non_cen_tot = 4, 4 annular = 5, 5 total = 6) c..ratio of earth's polar to equatorial radius double precision fac parameter (fac = 0.996633d0) c..number of earth equatorial radii in 1 au double precision auearth parameter (auearth = 149597870.0d0/6378.14d0) c..equatorial diameter of moon in earth equatorial radii double precision diam parameter (diam = 0.5450d0) c..equatorial diameter of sun in earth equatorial radii double precision dias parameter (dias = 218.25d0) c..popular format statements 222 format(1x,1p6e14.6) c..converting distances to earth radii and stretch z coordinates xm = xm1*auearth ym = ym1*auearth zm = zm1*auearth/fac rm = rm1*auearth xs = xs1*auearth ys = ys1*auearth zs = zs1*auearth/fac rs = rs1*auearth c..sun to moon unit vector xms = xm - xs yms = ym - ys zms = zm - zs dms = sqrt(xms**2 + yms**2 + zms**2) ex = xms/dms ey = yms/dms ez = zms/dms c..distance from moon to fundamental plane c..distance from earth to shadow axis s0 = -(xm*ex + ym*ey + zm*ez) delta = s0**2 + 1.0d0 - xm**2 - ym**2 - zm**2 r0 = sqrt(1.0d0 - delta) c..diameter of umbra and penumbra on the fundamental plane dumbra = (dias - diam) * s0/dms - diam dpenumbra = (dias + diam) * s0/dms + diam c..shadow axis intersects earth, so an annular or total eclipse if (r0 .lt. 1.0) then s = s0 - sqrt(delta) dumbra = (dias - diam) * s/dms - diam xshad = xm + s*ex yshad = ym + s*ey zshad = zm + s*ez zshad = zshad * fac if (dumbra .gt. 0.0) then iphase = annular else iphase = total end if c..non-central eclipses else if ( r0 .lt. (1.0d0 + 0.5d0*abs(dumbra)) ) then if (dumbra .gt. 0) then iphase = non_cen_ann else iphase = non_cen_tot end if c..partial eclipse else if ( r0 .lt. (1.0d0 + 0.5d0*dpenumbra) ) then iphase = partial c..no eclipse else iphase = noeclipse end if return end subroutine sunmoon(tjd,xm,ym,zm,rm,xs,ys,zs,rs) implicit none save c..simply returns the moon and suns equatorial postion vectors c..declare double precision tjd,xm,ym,zm,rm,xs,ys,zs,rs c..for getting the positional data character*40 aname integer isun,imoon,icoord parameter (isun=1, imoon=2, icoord = 1) double precision glon,glat,height, 1 xequ,yequ,zequ,vxequ,vyequ,vzequ, 2 ra,dec,rah,ram,ras,decd,decm,decs, 3 xecl,yecl,zecl,vxecl,vyecl,vzecl, 4 lam,bet,lamd,lamm,lams,betd,betm,bets, 5 xalt,yalt,zalt,vxalt,vyalt,vzalt, 6 alt,azi,altd,altm,alts,azid,azim,azis, 7 dist,drdt,dsun,dsundt, 8 xmanom,semi_maj,qper,xecen,lan,xinc,aop, 9 period,tjdp parameter (glon = 0.0d0, 1 glat = 0.0d0, 2 height = 0.0d0) c..get the position of the sun call position( 1 isun,tjd,icoord,glon,glat,height, 2 aname, 3 xequ,yequ,zequ,vxequ,vyequ,vzequ, 4 ra,dec,rah,ram,ras,decd,decm,decs, 5 xecl,yecl,zecl,vxecl,vyecl,vzecl, 6 lam,bet,lamd,lamm,lams,betd,betm,bets, 7 xalt,yalt,zalt,vxalt,vyalt,vzalt, 8 alt,azi,altd,altm,alts,azid,azim,azis, 9 dist,drdt,dsun,dsundt, & semi_maj,qper,xecen,lan,xinc,aop,xmanom,period,tjdp) xs = xequ ys = yequ zs = zequ rs = dist c..get the position of the moon call position( 1 imoon,tjd,icoord,glon,glat,height, 2 aname, 3 xequ,yequ,zequ,vxequ,vyequ,vzequ, 4 ra,dec,rah,ram,ras,decd,decm,decs, 5 xecl,yecl,zecl,vxecl,vyecl,vzecl, 6 lam,bet,lamd,lamm,lams,betd,betm,bets, 7 xalt,yalt,zalt,vxalt,vyalt,vzalt, 8 alt,azi,altd,altm,alts,azid,azim,azis, 9 dist,drdt,dsun,dsundt, & semi_maj,qper,xecen,lan,xinc,aop,xmanom,period,tjdp) xm = xequ ym = yequ zm = zequ rm = dist return end include 'position_routine.f' include 'glue.f'