procedure mkdifflc # Extract light curve of a SN from an image stack using image subtraction

int    beginste  {1,min=1,max=6,prompt="The beginning step(1=shift,2=register,3=diff,4=errcalc)"}
string imagelist {'imagelist', prompt="list of images"}
int    starnum   {10,prompt="number of stars for registration"}
real   ximsize   {2048, prompt="image size in pixels, x"}
real   yimsize   {2048, prompt="image size in pixels, y"}
real   psfsize   {21, prompt="psf size to use in fcpmxn"}
int    diffnum   {prompt="number of independent image (excluding the reference)"}

struct *lis
struct *temp
struct *temp1

begin

int i,j
int bgin
int num
string imname
real x1,y1
int xmin,ymin,xmax,ymax
real xref,yref,xnew,ynew
string junk
string refimname
string yesno,yesno1
real psfwidth
int delta
int arty1,arty2,arty3,arty4,arty5,arty6
int refmag
real refsky
real xrefstar,yrefstar
real xsn,ysn
int errorstarnum
real xfake[10],yfake[10]
real dmag[9],dflux[9]
real magsn,fluxsn,magrefstar,fluxrefstar
real magfake
real dfakemag[10],dfakeflux[10]
real magsum, magoffset
real tempmag,tempflux
string imdate

bgin = beginste

if (bgin == 1){
   goto stageshift
}
if (bgin == 2){
   goto stageregister
}
if (bgin == 3){
   goto stagediff
}
if (bgin == 4){
   goto stageerrorcalc
}

goto iend


stageshift: # rough alignment of the images

delete('shifts')

lis=imagelist

num=1

print ('mark object with comma, then q in each image')

while (fscan (lis, imname) != EOF){

display (imname,1)

delete('junk')

imexam (input=imname,frame=1, >> 'junk')

temp='junk'

j = fscan(temp,junk)
j = fscan(temp,x1,y1)

delete('junk')

if (num == 1){

   xref=x1
   yref=y1

   print (0,' ',0, >> 'shifts')

}

if (num > 1){
 
   xnew=xref-x1
   ynew=yref-y1

   print (xnew,' ',ynew, >> 'shifts')

}

num = num + 1

}

print('shifting')

lis=imagelist
temp='shifts'

while (fscan (lis, imname) != EOF){

      i = fscan (temp,x1,y1)

      delete(imname//'.shift.fits')

      imshift(input=imname,output=imname//'.shift',xshift=x1,yshift=y1)

}

stageregister: # register the images using geomap and geotran

print ('register')

lis=imagelist

num=1

while (fscan (lis, imname) != EOF){

      if (num == 1){ 

         delete ('geolist')

         display (imname//'.shift',1)

         print ('select stars')

         i = 1

         while (i<starnum+1){

               delete('junk')

               imexam (input=imname//'.shift',frame=1, >> 'junk')

               temp='junk'

               j = fscan(temp,junk)
               j = fscan(temp,x1,y1)

               print(i)
               print (x1,y1,i,>>'geolist')

               delete ('temptvcoo')
               print (x1,y1,>>'temptvcoo')
               tvmark (frame=1,coords='temptvcoo',mark='circle',radii='5,10,15,20',color=205)

               i=i+1

         }

         i=i-1

         print ('done: ',i,' stars selected.')

      }

      # for num > 1

      if (num > 1) {

      delete (imname//'.geo')

      temp='geolist'

      while (fscan (temp, xref, yref) != EOF){

            delete('tempcoo')
            delete('text')
            print(xref,' ',yref, >> 'tempcoo')

            imexamine(input=imname//'.shift',frame=1,image=imname//'.shift',logfile="imexa",defkey=",", imagecu="tempcoo",use_dis=no, >> 'text')
            temp1='text'
            i = fscan(temp1, junk)
            i = fscan(temp1, xnew,ynew)
            print(xref,' ',yref,' ',xnew,' ',ynew, >> imname//'.geo')

      }

      }

      num=num+1

}

lis=imagelist

num=1

while (fscan (lis, imname) != EOF){

      if (num == 1){

         delete(imname//'.reg.fits')

         imcopy (imname//'.shift',imname//'.reg')

      }

      if (num > 1){

         delete(imname//'.reg.fits')

         geomap(input=imname//'.geo',database="database",xmin=1,xmax=ximsize,ymin=1,ymax=yimsize,transform=imname,fitgeom="general",function="polynomial",xxorder=2,xyorder=2,yyorder=2,yxorder=2,interact=yes)

         geotran(input=imname//'.shift',output=imname//'.reg',database="database",transform=imname,geometr="geometric",xmin=1,xmax=ximsize,ymin=1,ymax=yimsize)

      }

      num = num + 1

}

stagediff: # difference all the images against a reference

# trim images (optional)

print ('do you want to trim the images before differencing ? (y/n)')

scan yesno

if (yesno == 'y'){

   lis=imagelist
   i = fscan (lis, imname)

   remark:

   display (imname,1)
   print ('enter limits: xmin')
   scan (xmin)
   print ('enter limits: xmax')
   scan (xmax)
   print ('enter limits: ymin')
   scan (ymin)
   print ('enter limits: ymax')
   scan (ymax)

   delete('boxcoo')

   print(xmin,' ',ymin,>> 'boxcoo')
   print(xmin,' ',ymax,>> 'boxcoo')
   print(xmin,' ',ymax,>> 'boxcoo')
   print(xmax,' ',ymax,>> 'boxcoo')
   print(xmax,' ',ymax,>> 'boxcoo')
   print(xmax,' ',ymin,>> 'boxcoo')
   print(xmax,' ',ymin,>> 'boxcoo')
   print(xmin,' ',ymin,>> 'boxcoo')

   tvmark(frame=1,coords='boxcoo',mark='line')

   print ('is this ok? (y/n)')

   scan yesno1

   if (yesno1 == 'n'){goto remark}

   ximsize=xmax-xmin+1
   yimsize=ymax-ymin+1

   print ('trimming')

   lis=imagelist

   while (fscan (lis, imname) != EOF){ 
  
         imdel(imname//'.reg.big.fits')
         imrename(imname//'.reg.fits',imname//'.reg.big.fits')
         imcopy(imname//'.reg.big.fits['//xmin//':'//xmax//','//ymin//':'//ymax//']', imname//'.reg.fits')

   }

}

# in case images already trimmed by previous runs

if (yesno == 'n'){

   print ('is ximsize='//ximsize//' ? (y/n)')
   scan (yesno1)
   if (yesno1 == 'n'){
      print ('enter new ximsize')
      scan ximsize
   }

   print ('is yimsize='//yimsize//' ? (y/n)')
   scan (yesno1)
   if (yesno1 == 'n'){
      print ('enter new yimsize')
      scan yimsize
   }

}

lis=imagelist

delete ('difflist')

num = 1

while (fscan (lis, imname) != EOF){

      if (num == 1){
 
         # selecting PSF stars

         print ('image ',imname//'.reg',' selected as reference.')

         refimname=imname//'.reg'

         display (refimname,1)

         print ('select PSF stars')
         print ('to end selection, hit same star twice')

         i = 1

         # these are used to track the last star and stop the selection

         xref=0
         yref=0

         delete ('PSFstarlist')

         while (i<20){ # arbitrary maximum for number of PSF stars

               delete('junk')

               imexam (input=refimname,frame=1, >> 'junk')

               temp='junk'

               j = fscan(temp,junk)
               j = fscan(temp,x1,y1,junk,junk,junk,junk,junk,junk,junk,junk,psfwidth)

               # if same star selected twice - exit the loop

               if (x1==xref){
                  if (y1==yref){
                     goto PSFend
                  }
               }

               print(i)
               print (x1,y1,>>'PSFstarlist')

               xref=x1
               yref=y1

               delete ('temptvcoo')
               print (x1,y1,>>'temptvcoo')
               tvmark (frame=1,coords='temptvcoo',mark='circle',radii='5,10,15,20',color=205)

               i=i+1
         }

         PSFend:

         i=i-1

         print ('done: ',i,' stars selected.')

         # set flux and insert reference star

         delta=abs(ximsize/8)
         delete ('artstars')
         arty1=yimsize-delta-1
         arty2=yimsize-2*delta-1
         arty3=yimsize-3*delta-1
         arty4=yimsize-4*delta-1
         arty5=yimsize-5*delta-1
         arty6=yimsize-6*delta-1

         print (delta,' ',arty1,' ',10, >> 'artstars')
         print (delta,' ',arty2,' ',12, >> 'artstars')
         print (delta,' ',arty3,' ',14, >> 'artstars')
         print (delta,' ',arty4,' ',16, >> 'artstars')
         print (delta,' ',arty5,' ',18, >> 'artstars')
         print (delta,' ',arty6,' ',20, >> 'artstars')

         delete('reftest.fits')
         psfwidth=psfwidth/2 # use the seeing of the last psf star to set the artificial seeing
         mkobjects(input=refimname,output='reftest',objects='artstars',magzero=25,star='gaussian',radius=psfwidth)
         display('reftest')
 
         print('select bright, unsaturated artificial source')
 
         delete('junk')
         imexam (input='reftest',frame=1, >> 'junk')
         temp='junk'
         j = fscan(temp,junk)
         j = fscan(temp,x1,y1)

         if (abs(y1-arty1) < delta/2) {refmag=10}
         if (abs(y1-arty2) < delta/2) {refmag=12}
         if (abs(y1-arty3) < delta/2) {refmag=14}
         if (abs(y1-arty4) < delta/2) {refmag=16}
         if (abs(y1-arty5) < delta/2) {refmag=18}
         if (abs(y1-arty6) < delta/2) {refmag=20}

         # set area to determine sky value in reference image (optional)

         print ('would you like to set the sky area ? (y/n)')
         
         scan (yesno)
 
         if (yesno == 'y'){

            print ('set xmin:')
            scan (xmin)
            print ('set xmax:')
            scan (xmax)
            print ('set ymin:')
            scan (ymin)
            print ('set ymax:')
            scan (ymax)

         }

         if (yesno == 'n'){
            xmin=1
            xmax=ximsize
            ymin=1
            ymax=yimsize

         }

         delete ('junk')
         imstat(images=refimname//'['//xmin//':'//xmax//','//ymin//':'//ymax//']',fields='midpt', >> 'junk')
         temp='junk'
         j = fscan(temp,junk)
         j = fscan(temp,refsky)

         # inserting artificial reference star into reference image

         imdel (refimname//'.nosky')
         imarith (refimname,'-',refsky,refimname//'.nosky')
         imdel(refimname//'.inv')
         imarith (refimname//'.nosky','*',-1,refimname//'.inv')
         display (refimname,1)
         print ('mark position for reference star, use x key')
         delete ('junk')
         imexam (input=refimname,frame=1, >> 'junk')
         temp='junk'
         j = fscan(temp,xrefstar,yrefstar)
         delete ('artrefstar')
         print (xrefstar,' ',yrefstar,' ',refmag, >> 'artrefstar')
         imdel(refimname//'.inv.art')
         mkobjects(input=refimname//'.inv',output=refimname//'.inv.art',objects='artrefstar',magzero=25,star='gaussian',radius=psfwidth) 
         imdel(refimname//'.art')
         imarith (refimname//'.inv.art','*',-1,refimname//'.art')
         imarith (refimname//'.art','+',refsky,refimname//'.art')
      }

      if (num > 1){

         # subtract reference from all images 

         fcpmxn(newimag=imname//'.reg',refimag=refimname//'.art',region='PSFstarlist',ximsize=ximsize,yimsize=yimsize,psfsize=psfsize,statsiz=50,datamin=-5000,datamax=500000)
         print(imname//'.reg.dif.fits', >> 'difflist')

      }

      num=num+1

}

stageerrorcalc: # calculate error estimate for SN using artificial stars

lis='difflist'

# determining the SN coordinates

while (fscan (lis, imname) != EOF){ 
  
      display (imname,1)
      print ('is this subtractin good enough to mark SN ? (y/n)')
      scan (yesno)

      if (yesno == 'y'){

         print ('mark SN, use comma')

         delete ('junk')
         imexam (input=imname,frame=1, >> 'junk')
         temp='junk'
         j = fscan(temp,junk)
         j = fscan(temp,xsn,ysn)

         delete ('temptvcoo')
         print (xsn,ysn,>>'temptvcoo')
         tvmark (frame=1,coords='temptvcoo',mark='circle',radii='5,10,15,20',color=205)

         goto endsnselect
      
      }

}

endsnselect:

lis=imagelist

# setting fake star locations

j = fscan (lis, imname) 
  
display (imname//'.reg',1)

tvmark (frame=1,coords='temptvcoo',mark='circle',radii='5,10,15,20',color=206)

print ('how many fake stars to use for error estimation (up to 10) ?')

scan (errorstarnum)

print ('mark up to '//errorstarnum//' positions for fake SNe, use the x key')

i=0

while (i<errorstarnum){

       delete ('junk')
       imexam (input=imname,frame=1, >> 'junk')
       temp='junk'
       j = fscan(temp,x1,y1)

       xfake[i+1]=x1
       yfake[i+1]=y1

       delete ('temptvcoo')
       print (x1,y1,>>'temptvcoo')
       tvmark (frame=1,coords='temptvcoo',mark='circle',radii='5,10,15,20',color=205)

       print(i)

       i=i+1

}

# read position and magnitude of reference star

temp='artrefstar'
j = fscan(temp,xrefstar,yrefstar,refmag)

# go over difflist, measure SN-ref difference

lis=imagelist

delete ('errorcalclist')

num=1

delete ('tempcoo')
print (xsn,' ',ysn,>>'tempcoo')
print (xrefstar,' ',yrefstar,>>'tempcoo')

i=1

while (fscan (lis, imname) != EOF){ 

      if (num == 1){

         refimname=imname//'.reg' # defining refimname to be used in "iterate" stage

         goto skipreference # skip reference image

      } 

      # measure PSF in diff image from reference artificial star

      delete ('junk')
      imexamine(input=imname//'.reg.dif.fits',frame=1,image=imname//'.reg.dif.fits',logfile="imexa",defkey=",", imagecu="artrefstar",use_dis=no, >> 'junk')
      temp='junk'
      j = fscan(temp,junk)
      j = fscan(temp,junk,junk,junk,junk,junk,junk,junk,junk,junk,junk,psfwidth)

      rimexam.radius=psfwidth

      delete ('junk')
      imexamine(input=imname//'.reg.dif.fits',frame=1,image=imname//'.reg.dif.fits',logfile="imexa",defkey=",", imagecu="tempcoo",use_dis=no, >> 'junk')
      temp='junk'
      j = fscan(temp,junk)
      j = fscan(temp,junk,junk,magsn,fluxsn)
      j = fscan(temp,junk,junk,magrefstar,fluxrefstar)
      dmag[i]=magsn-magrefstar
      dflux[i]=fluxsn/fluxrefstar

      # insert artificial reference stars

      j=1

      delete('fakestarlist')

      while (j<errorstarnum+1){

            magfake=refmag # as first guess to 
            print (xfake[j],' ',yfake[j],' ',magfake,>> 'fakestarlist')
            j=j+1
    
      } 

      imdel(imname//'.errorcalc')

      # get PSF in the source image from PSF stars 

      delete ('junk')
      imexamine(input=imname//'.reg.fits',frame=1,image=imname//'.reg.fits',logfile="imexa",defkey=",", imagecu="PSFstarlist",use_dis=no, >> 'junk')
      temp='junk'
      j = fscan(temp,junk)
      j = fscan(temp,junk,junk,junk,junk,junk,junk,junk,junk,junk,junk,psfwidth)

      rimexam.radius=psfwidth

      # divide by 2 to get radius of PSF

      psfwidth=psfwidth/2 # change FWHM to radius

      mkobjects(input=imname//'.reg.fits',output=imname//'.errorcalc',objects='fakestarlist',magzero=25,star='gaussian',radius=psfwidth) 

      print (imname//'.errorcalc', >> 'errorcalclist')

      i=i+1

      skipreference:

      num=num+1

}

# iterate

delete ('tempcoo')
print (xsn,' ',ysn,>>'tempcoo')
print (xrefstar,' ',yrefstar,>>'tempcoo')

i=1

while (i<errorstarnum+1){
      print (xfake[i],' ',yfake[i],>> 'tempcoo')
      i=i+1
}

lis='errorcalclist'
delete('errordifflist')


while (fscan (lis, imname) != EOF){

      # subtract reference from all errorcalc images 

      fcpmxn(newimag=imname,refimag=refimname//'.art',region='PSFstarlist',ximsize=ximsize,yimsize=yimsize,psfsize=psfsize,statsiz=50,datamin=-5000,datamax=500000)

}

# measure offest between SN and fake stars

lis=imagelist

delete ('lclist')

num=1

while (fscan (lis, imname) != EOF){ 

      if (num == 1){

         refimname=imname//'.reg' # defining refimname to be used below

         goto skipreference2 # skip reference image

      } 

      # set measurement radius in diff image by the artificial reference star

      delete ('junk')
      imexamine(input=imname//'.errorcalc.dif.fits',frame=1,image=imname//'.errorcalc.dif.fits',logfile="imexa",defkey=",", imagecu="artrefstar",use_dis=no, >> 'junk')
      temp='junk'
      j = fscan(temp,junk)
      j = fscan(temp,junk,junk,junk,junk,junk,junk,junk,junk,junk,junk,psfwidth)

      rimexam.radius=psfwidth

      # measure offset between mean of fake stars and SN

      delete ('junk')
      imexamine(input=imname//'.errorcalc.dif.fits',frame=1,image=imname//'.errorcalc.dif.fits',logfile="imexa",defkey=",", imagecu="tempcoo",use_dis=no, >> 'junk')
      temp='junk'
      j = fscan(temp,junk)
      j = fscan(temp,junk,junk,magsn,fluxsn)
      j = fscan(temp,junk) # skip line with artificial reference star which is not useful here
      i=1
      while (i<errorstarnum+1){
            j = fscan(temp,junk,junk,magrefstar,fluxrefstar)
            dfakemag[i]=magsn-magrefstar
            dfakeflux[i]=fluxsn/fluxrefstar
            i=i+1
      }

      # calculate mean mag offset

      i=1
      magsum=0
      while (i<errorstarnum+1){
            magsum=magsum+dfakemag[i]
            i=i+1
      }

      magoffset=magsum/errorstarnum

      j=1

      delete('fakestarlist')

      while (j<errorstarnum+1){

            magfake=refmag+magoffset # now accurate fake mag
            print (xfake[j],' ',yfake[j],' ',magfake,>> 'fakestarlist')
            j=j+1
    
      } 

      imdel(imname//'.errorcalc.2')

      # measure PSF in source image using PSF stars

      delete ('junk')
      imexamine(input=imname//'.reg.fits',frame=1,image=imname//'.reg.fits',logfile="imexa",defkey=",", imagecu="PSFstarlist",use_dis=no, >> 'junk')
      temp='junk'
      j = fscan(temp,junk)
      j = fscan(temp,junk,junk,junk,junk,junk,junk,junk,junk,junk,junk,psfwidth)

      rimexam.radius=psfwidth

      # divide by 2 to get radius of PSF

      psfwidth=psfwidth/2 # change FWHM to radius

      # second iteration simulation and subtraction

      mkobjects(input=imname//'.reg.fits',output=imname//'.errorcalc.2',objects='fakestarlist',magzero=25,star='gaussian',radius=psfwidth) 

      fcpmxn(newimag=imname//'.errorcalc.2.fits',refimag=refimname//'.art',region='PSFstarlist',ximsize=ximsize,yimsize=yimsize,psfsize=psfsize,statsiz=50,datamin=-5000,datamax=500000)

      # final measuring of the SN, reference star and fake stars

      delete (imname//'.finallc')

      # set measurement radius in diff image by artificial star

      delete ('junk')
      imexamine(input=imname//'.errorcalc.2.fits.dif.fits',frame=1,image=imname//'.errorcalc.2.fits.dif.fits',logfile="imexa",defkey=",", imagecu="artrefstar",use_dis=no, >> 'junk')
      temp='junk'
      j = fscan(temp,junk)
      j = fscan(temp,junk,junk,junk,junk,junk,junk,junk,junk,junk,junk,psfwidth)

      rimexam.radius=psfwidth

      # get date from header

      imgets (image=imname,param='obsmjd')

      imdate=imgets.value # for old ccd

      if (imgets.value == '0'){

         # newccd has different keyword

         imgets (image=imname,param='obsjd')

         imdate=imgets.value

      }

      # measure offset between mean of fake stars and SN

      delete ('junk')
      imexamine(input=imname//'.errorcalc.2.fits.dif.fits',frame=1,image=imname//'.errorcalc.2.fits.dif.fits',logfile="imexa",defkey=",", imagecu="tempcoo",use_dis=no, >> 'junk')
      temp='junk'
      j = fscan(temp,junk)
      j = fscan(temp,junk,junk,magsn,fluxsn)
      j = fscan(temp,junk,junk,magrefstar,fluxrefstar)

      tempmag=magsn-magrefstar
      print ('% magnitude offset from the reference artificial star, and JD/MJD:',>>imname//'.finallc')
      print(tempmag,' ',imdate, >>imname//'.finallc')
      
      tempflux=fluxsn/fluxrefstar
      print ('% ratio of SN flux to artificial star flux, and JD/MJD:',>>imname//'.finallc')
      print(tempflux,' ',imdate,>>imname//'.finallc')

      print ('% SN magnitude and flux:',>>imname//'.finallc')
      print (magsn,' ',fluxsn,>>imname//'.finallc')
   
      print ('% fake stars magnitudes and fluxes:',>>imname//'.finallc')

      i=1
      while (i<errorstarnum+1){
            j = fscan(temp,junk,junk,magrefstar,fluxrefstar)
            print (magrefstar,' ',fluxrefstar,>>imname//'.finallc')
            i=i+1
      }

      print (imname//'.finallc',>> 'lclist')

      skipreference2:

      num=num+1

}


iend: # end of script
print("Ending...  Bye.")
beep

end