*23456789a123456789b123456789c123456789d123456789e123456789f123456789g12 * program lc_fit parameter (npp=3,np=150,nxx=3,nx=1000) parameter (nyy=np*nxx,nyyy=nyy*nxx) dimension ip(npp),timep(np,npp),xmagp(np,npp), & ix(nxx),timex(nx,nxx),xmagx(nx,nxx), & xoffset(nxx),xmago(nxx),time(np) dimension iorder(nyy),iinteger(nyy),timexx(nyy),xmagxx(nyy,nxx) character cchar(nyy)*1 character filenp(npp)*30,filenx(nxx)*30 data filenp/'./lc_brunou.dat','./lc_brunob.dat', & './lc_brunov.dat'/ * data filenx/'./sn1972e/lc_u.dat','./sn1972e/lc_b.dat', * & './sn1972e/lc_v.dat'/ * data filenx/'./sn1984a/lc_u.dat','./sn1984a/lc_b.dat', * & './sn1984a/lc_v.dat'/ * data filenx/'./sn1986g/lc_u.dat','./sn1986g/lc_b.dat', * & './sn1986g/lc_v.dat'/ * data filenx/'./sn1981b/lc_u.dat','./sn1981b/lc_b.dat', * & './sn1981b/lc_v.dat'/ * data filenx/'./sn1990n/lc_u.dat','./sn1990n/lc_b.dat', * & './sn1990n/lc_v.dat'/ * data filenx/'./sn1992a/lc_u.dat','./sn1992a/lc_b.dat', * & './sn1992a/lc_v.dat'/ data filenx/'./sn1994d/lc_u.dat','./sn1994d/lc_b.dat', & './sn1994d/lc_v.dat'/ * data filenx/'./sn1994ae/lc_u.dat','./sn1994ae/lc_b.dat', * & './sn1994ae/lc_v.dat'/ * data filenx/'./sn1995d/lc_u.dat','./sn1995d/lc_b.dat', * & './sn1995d/lc_v.dat'/ data ix/nxx*0/,xmagxx/nyyy*0./ * * The isuper is the year*1000 plus the order number of the letter * designation. * a b c d e f g h i j k l m n * 1 2 3 4 5 6 7 8 9 10 11 12 13 14 * * SN 1972E isuper=1972000+5 * isuper=1972000+5 * SN 1981B isuper=1981000+2 * isuper=1981000+2 * SN 1984A isuper=1984000+1 * isuper=1984000+1 * SN 1986G isuper=1984000+1 * isuper=1986000+7 * SN 1990N isuper=1990000+14 * isuper=1990000+14 * SN 1992a isuper=1992000+1 * isuper=1992000+1 * SN 1994d isuper=1994000+4 isuper=1994000+4 * SN 1994ae isuper=1994000+31 * isuper=1994000+31 * SN 1995d isuper=1995000+4 * isuper=1995000+4 * * Read in Bruno's template light curves. * do 410 i=1,npp open(unit=1,file=filenp(i),status='old') * call readhead(1,0,0) * j=1 120 continue read(1,*,end=130) timep(j,i),xmagp(j,i) if(isuper .eq. 1990000+14 .and. i .eq. 3) then * This is just because the V maximum of 90N occurred only * one day later than the B maximum. So we shift the V * template to 1 day earlier. timep(j,i)=timep(j,i)-1 end if * print*,j,timep(j,i),xmagp(j,i) j=j+1 go to 120 130 continue ip(i)=j-1 close(unit=1) * 410 continue * * * Read in the observed light curves. * k=0 do 420 i=1,nxx open(unit=2,file=filenx(i),status='unknown') * j=1 read(2,*,err=160,end=160) call readhead(2,0,0) 150 continue read(2,*,err=160,end=160) timex(j,i),xmagx(j,i) * print*,timex(j,i),xmagx(j,i) if(isuper .eq. 1992000+1) then timex(j,i)=timex(j,i)-48622.7 else if(isuper .eq. 1994000+4) then timex(j,i)=timex(j,i)-49432.90+18. * This is Vacca & Leibundgut 1996 B max. else if(isuper .eq. 1994000+31) then timex(j,i)=timex(j,i)-9666.5 else if(isuper .eq. 1995000+4) then timex(j,i)=timex(j,i)-9751.5 end if k=k+1 timexx(k)=timex(j,i) xmagxx(k,i)=xmagx(j,i) * print920,timex(j,i),xmagx(j,i) * print920,timex(j,i),0.,0.,xmagx(j,i),0.,0. j=j+1 go to 150 160 continue ix(i)=j-1 close(unit=2) * 420 continue * itimexx=k * call ordering(2,itimexx,iorder,iinteger,timexx,cchar) do 421 i=1,nxx call ordering(-2,itimexx,iorder,iinteger,xmagxx(1,i),cchar) 421 continue i=0 do 422 k=1,itimexx * print920,k,timexx(k),(xmagxx(k,m),m=1,nxx) if(timexx(k) .gt. 0.) then i=i+1 timexx(i)=timexx(k) do 423 m=1,nxx xmagxx(i,m)=xmagxx(k,m) 423 continue do 424 l=k+1,itimexx if(timexx(l) .gt. timexx(i)+1.e-5) go to 170 timexx(l)=-1. do 426 m=1,nxx if(abs(xmagxx(l,m)) .gt. 1.e-5 .and. & abs(xmagxx(i,m)) .lt. 1.e-5 ) then xmagxx(i,m)=xmagxx(l,m) else if(abs(xmagxx(l,m)) .gt. 1.e-5) then xmagxx(i,m)=.5*(xmagxx(i,m)+xmagxx(l,m)) * I'm assuming there will never be any more than two * separate observations with numerically identical * times. If there are, the error is small: I'm just * doing an only slightly stupid unequal averaging. end if 426 continue 424 continue 170 continue end if 422 continue itimexx=i * print* print908 908 format('Photometry: observed') do 428 j=1,itimexx print920,timexx(j),(xmagxx(j,i),i=1,nxx),0.,0. 428 continue print930 * * Now find the mean offset between the templates and the * observed light curves. * print* do 430 i=1,nxx if(ix(i) .eq. 0) go to 430 sum=0. icount=0 do 440 j=1,ix(i) * if(isuper .eq. 1981000+2 .and. ! This is for * & (timex(j,i) .lt. 14. .or. ! a maximum light fit. * & timex(j,i) .gt. 22. )) go to 440 if(isuper .eq. 1994000+4 .and. ! This is for & (timex(j,i) .lt. 40. .or. ! an early decline fit. & timex(j,i) .gt. 100. )) go to 440 icount=icount+1 call interpo(1,timex(j,i),timep(1,i),xmagp(1,i), & ip(i),1,1,xmagpp) * print*,i,timex(j,i),xmagx(j,i),xmagpp sum=sum+(xmagx(j,i)-xmagpp) 440 continue xoffset(i)=sum/real(icount) print*,'The ',i,'th offset magnitude is' print*,xoffset(i) print* 430 continue * * Read in the full set of Bruno's times two late Branch & Doggett * times. * open(unit=3,file='./lc_brunot.dat',status='old') read(3,*) itime do 450 i=1,itime read(3,*) time(i) 450 continue * * Print out the fitted light curve. * print910 910 format('Photometry: fitted') do 460 i=1,itime do 470 j=1,nxx if(ix(j) .gt. 0) then call interpo(1,time(i),timep(1,j),xmagp(1,j), & ip(j),1,1,xmagpp) xmago(j)=xmagpp+xoffset(j) else xmago(j)=0. end if 470 continue print920,time(i),(xmago(j),j=1,nxx),0.,0. * print920,time(i),xmago(3) 460 continue * 920 format(f10.3,5f10.4) print930 930 format('END:') * delx=1. do 480 i=1,nxx k=0 do 482 j=1,itimexx if(abs(xmagxx(j,i)) .gt. 1.e-5) then k=k+1 timex(k,i)=timexx(j) xmagx(k,i)=xmagxx(j,i) iorder(k)=j end if 482 continue itimex=k * * smooth over 1 day. Remember this means half DELX ahead * and half DELX behind. call smooth(itimex,timex(1,i),xmagx(1,i),delx,0.) do 484 j=1,itimex xmagxx(iorder(j),i)=xmagx(j,i) print920,timex(j,i),xmagx(j,i) 484 continue * 480 continue * print* print940,delx 940 format('Photometry: observed: box-car smoothed', & ' over a ',f5.2,' day interval') do 500 j=1,itimexx,1 print920,timexx(j),(xmagxx(j,i),i=1,nxx),0.,0. 500 continue print930 * end