soimosaic.cl

procedure soimosaic (inimage)

# Mosaic the  extensions of the SOAR Optical Imager

char inimage	{"",prompt="Input image or list"}
char outimage   {"",prompt="Output image or list"}
char outpref    {"m",prompt="Prefix for output images"}     
real xgap	    {102.,prompt="Size of the gap between the CCD's (unbinned pixels)"}
real ygap	    {6.,prompt="Difference in y between the CCD's (unbinned pixels)"}
real tilt       {0.,prompt="Misalignment angle between the CCD's (degrees)"}
bool clobber    {no,prompt="Modify existing output images?"}


begin

bool merged,imcopyveb, l_clobber
char l_inimage,l_outimage,l_outpref
char outim,tmpinlist,tmpoutlist, tmpim[6],int_image,datsec,aux
char imname[1000], outname[1000]
int nc, nl, nim,nout,n_ext
real l_xgap, l_ygap, l_tilt, Xbin,Ybinr, rot[4]


# Verify if all needed packages has been loaded

if (deftask("proto")==no) {
  print("ERROR: The PROTO package is required but not defined.")
  goto crash
}

if (deftask("images")==no) {
  print("ERROR: The PROTO package is required but not defined.")
  goto crash
}

imcopyveb = imcopy.verbose
imcopy.verbose = no


#Expand the input image list

l_inimage = inimage
tmpinlist=mktemp("tmpinlist")
sections(l_inimage,opt="fullname", >tmpinlist)
nim = sections.nimages
if(nim>1000) {
  print("Number of images has to be smaller that 1000")
  goto crash
}

# Verify the existence of the input files

list = tmpinlist
for(i=1; fscan(list,s1) !=EOF; i+=1) {
  imname[i] = s1
  j = strlen(imname[i])
  s2 = substr(imname[i],j-4,j)
  if(s2 != ".fits") {
    imname[i]=imname[i]//".fits"
  }  
  if(!access(imname[i])) {
      printf ("\nERROR: Input image %s not found\n",imname[i])
	  goto crash
  } 
}
delete(tmpinlist,ver-)

#Expand the output image list

l_outimage = outimage
l_outpref = outpref

tmpoutlist=mktemp("tmpoutlist")

if (l_outimage=="") {
  #prefix
  for(i=1;i<=nim;i+=1)
	outname[i] = l_outpref//imname[i]
  }
else {
  sections(l_outimage,opt="fullname", >tmpoutlist)
  nout = sections.nimages
  if(nim != nout) {
    print("ERROR: Number of input and output images should be the same")
    goto crash
  } 	 
list = tmpoutlist
for(i=1; fscan(list,s1) !=EOF; i+=1) {
  outname[i] = s1
  j = strlen(outname[i])
  s2 = substr(outname[i],j-4,j)
  if(s2 != ".fits") 
	outname[i]=outname[i]//".fits"
  }
delete(tmpoutlist,ver-) 	
}  


for(k=1;k<=nim;k+=1) # Loop over the images
  {

  if(access(outname[k]) ) {
    if(clobber) {
      delete(outname[k])
    } else {
    printf ("\nERROR: Operation would overwrite existing image (%s) \n",outname[k])
	goto crash
    }	
  }

  #Create temporary files
  tmpim[1]=mktemp("tmpim1")
  tmpim[2]=mktemp("tmpim2")
  tmpim[3]=mktemp("tmpim3")
  tmpim[4]=mktemp("tmpim4")
  tmpim[5]=mktemp("tmpim5")
  

  # Test if amplifiers of the same CCD have been merged or not 
  imextensions (imname[k],lindex-,lname-,lver-, > "dev$null")
  n_ext = imextensions.nimages
  printf ("%s %d -> %s\n", imname[k], n_ext, outname[k])
  
  if(n_ext != 2 && n_ext != 4)
    {
	print("ERROR: Number of extensions should be 2 or 4 !")
	goto crash
	}

  if(n_ext == 2)
    merged = yes
  else
    merged = no

  if(!merged)
    {			
	hselect (imname[k]//"[1]", "naxis1", yes) | scan (nc)
	s1 = imname[k]//"[1],"//imname[k]//"[2]"
	imcombine (s1, tmpim[1], offsets="grid 2 "//nc,
	   headers="", bpmasks="", rejmasks="", nrejmasks="", expmasks="",
	   sigmas="", logfile="", combine="average", reject="none",
	   project=no, outtype="real", outlimits="", 
		masktype="none",
	   maskvalue=0, blank=0., scale="none", zero="none",
	   weight="none", statsec="", expname="", lthreshold=INDEF,
	   hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes,
	   lsigma=3., hsigma=3., rdnoise="0.", gain="1.", snoise="0.",
	   sigscale=0.1, pclip=-0.5, grow=0.)

	hselect (imname[k]//"[3]", "naxis1", yes) | scan (nc)
	s1 = imname[k]//"[3],"//imname[k]//"[4]"
	imcombine (s1, tmpim[2], offsets="grid 2 "//nc,
	   headers="", bpmasks="", rejmasks="", nrejmasks="", expmasks="",
	   sigmas="", logfile="", combine="average", reject="none",
	   project=no, outtype="real", outlimits="", 
		masktype="none",
	   maskvalue=0, blank=0., scale="none", zero="none",
	   weight="none", statsec="", expname="", lthreshold=INDEF,
	   hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes,
	   lsigma=3., hsigma=3., rdnoise="0.", gain="1.", snoise="0.",
	   sigscale=0.1, pclip=-0.5, grow=0.)
    }

  hselect(imname[k]//"[1]","CCDSUM",yes) | scan(s1,s2)

  if(strlen(s1) < 1)
    {
    printf ("\nERROR: Keyword CCDSUM not found on image %s\n",imname[i])
	goto crash
    }
		
  aux = substr(s1,2,2)
  Xbin = real (aux)
  aux = substr(s2,1,1)
  Ybin = real (aux)

  l_xgap = xgap / Xbin
  l_ygap = ygap / Ybin

  # Shifts in the y direction and rotates the right CCD

  geotran(tmpim[2],tmpim[4],"","",
    xshift=0.,yshift=l_ygap, xrotation= 0.,yrotation= 0., 
    xmag=1,ymag=1,xmin=INDEF,xmax=INDEF,ymin=INDEF,ymax=INDEF,
    ncols=INDEF,nlines=INDEF,fluxconserve=yes, 
    nxblock=2048,nyblock=2048,interpolant="linear",boundary="constant",
    constan=0,verbose=no)
  
  imcopy(tmpim[1],tmpim[3])
 
  hselect (tmpim[3], "naxis1", yes) | scan (nc); nc = nc + l_xgap
  s1 = tmpim[3]//","//tmpim[4]
  imcombine (s1, outname[k], offsets="grid 2 "//nc,
     headers="", bpmasks="", rejmasks="", nrejmasks="", expmasks="",
     sigmas="", logfile="", combine="average", reject="none",
	   project=no, outtype="real", outlimits="", 
		masktype="none",
	   maskvalue=0, blank=0., scale="none", zero="none",
     weight="none", statsec="", expname="", lthreshold=INDEF,
     hthreshold=INDEF, nlow=1, nhigh=1, nkeep=1, mclip=yes,
     lsigma=3., hsigma=3., rdnoise="0.", gain="1.", snoise="0.",
     sigscale=0.1, pclip=-0.5, grow=0.)
    
 # Delete useless keywords
   
  hedit(outname[k],"IMAGEID", "",add-,addonly-,delete+,verify-,show-,update+)
  hedit(outname[k],"CCDNAME", "",add-,addonly-,delete+,verify-,show-,update+)
  
  delete(tmpim[1]//".fits")
  delete(tmpim[2]//".fits")
  delete(tmpim[3]//".fits")
  delete(tmpim[4]//".fits")

}  

imcopy.verbose = imcopyveb

end