# Copyright(c) 2003-2005 Association of Universities for Research in Astronomy, Inc.

procedure mipsplit(inimages,outimages)

# This routine extracts the individual polarimetry images and 
# writes them to simple FITS files.
#
# Version:  October 18, 2005  KV wrote original script 
#
#           There is a header keyword POLANGLE for each extension that gives the 
#           values of the polarimetry plate angle in degrees.
#
char  inimages     {"",prompt="Input Michelle polarimetry image(s)"}	# OLDP-1-input-primary-single-prefix=s
char  outimages    {"",prompt="Output image(s) base name"}              # OLDP-1-output
char  rawpath      {"",prompt="Path for in raw images"}      		# OLDP-4
int status         {0,prompt="Exit status: (0=good, >0=bad)"}	        # OLDP-4
struct* scanfile   {"",prompt="Internal use only"}   			# OLDP-4

begin

    char   l_inputimages, l_outputimages, l_filename,l_rawpath
    char   in[100],out[100],header,exheader,instrument,l_temp
    char   tmpfile,filename
    int    i,j,k,l,m,n,nnods,nwp[4],jwp,nbad
    int    nimages,maximages,noutimages,l_extensions
    int    modeflag,inod

    tmpfile=mktemp("tmpin")
    l_inputimages=inimages
    l_outputimages=outimages
    l_rawpath=rawpath

    nimages=0
    maximages=100
    status=0
    nwp[1]=0
    nwp[2]=0
    nwp[3]=0
    nwp[4]=0
 
    #----------------------------------------------------------------------
    # Check the rawpath name for a final /
    if (substr(l_rawpath,(strlen(l_rawpath)),(strlen(l_rawpath))) != "/")
        l_rawpath=l_rawpath//"/"
    if (l_rawpath=="/" || l_rawpath==" ")
        l_rawpath=""

    #-----------------------------------------------------------------------
    # Load up arrays of input name lists
    # This version handles both *s, / and commas in l_inputimages

    # check that list file exists
    if (substr(l_inputimages,1,1)=="@") {
        l_temp=substr(l_inputimages,2,strlen(l_inputimages))
        if (!access(l_temp) && !access(l_rawpath//l_temp)) {
	    print("Error: Input file "//l_temp//" not found.")
	    status=1
	    goto clean
        }
    }

    # Count the number of in images
    # First, generate the file list if needed

    if (stridx("*",l_inputimages) > 0) {
        files(l_inputimages, > tmpfile)
        l_inputimages="@"//tmpfile
    }

    if (substr(l_inputimages,1,1)=="@")
        scanfile=substr(l_inputimages,2,strlen(l_inputimages))
    else {
        files(l_inputimages,sort-, > tmpfile)
        scanfile=tmpfile
    }

    i=0

    while (fscan(scanfile,l_filename) != EOF && i <= 100) {

	i=i+1

	if (substr(l_filename,strlen(l_filename)-4,strlen(l_filename)) == ".fits")
	    l_filename=substr(l_filename,1,strlen(l_filename)-5)

	if (!imaccess(l_filename) && !imaccess(l_rawpath//l_filename)) {
	    print("Error: Input image"//l_filename//" was not found.")
	    status=1
	    goto clean
	} else {
	    nimages=nimages+1
	    if (nimages > maximages) {
              print("Maximum number of input images exceeded")
		status=1 
		goto clean
	    }

	    if (l_rawpath=="" || l_rawpath==" ")
                in[nimages]=l_filename
	    else
                in[nimages]=l_rawpath//l_filename
	}
    }

    scanfile=""
    delete(tmpfile,ver-,>& "dev$null")

    if (nimages == 0) {
        print("Error:  No input images defined.")
	status=1
	goto clean
    }

    # Now, do the same counting for the out file

    nbad=0
    noutimages=0
    if ((l_outputimages != "") && (l_outputimages != " ")) {
	if (substr(l_outputimages,1,1) == "@")
	    scanfile=substr(l_outputimages,2,strlen(l_outputimages))
	else {
	    if (stridx("*",l_outputimages) > 0) {
		files(l_outputimages,sort-) | \
		    match(".hhd",stop+,print-,metach-, > tmpfile)
		scanfile=tmpfile
	    } else {
		files(l_outputimages,sort-, > tmpfile)
		scanfile=tmpfile
	    }
	}
  
	while (fscan(scanfile,l_filename) != EOF) {
	    noutimages=noutimages+1
	    if (noutimages > maximages) {
	        print("Error: Maximum number of output images exceeded.")
		status=1
		goto clean
	    }
            if (substr(l_filename,strlen(l_filename)-4,strlen(l_filename)) != ".fits") 
              l_filename=l_filename//".fits"
# test for one of the output files, exit if it is there....
            out[noutimages]=substr(l_filename,1,strlen(l_filename)-5)//"_wp1im1.fits"
 	    if (imaccess(out[noutimages])) {
	       print("Error: Output image "//out[noutimages]//" already exits.")
               nbad+=1
	    }
# The output file name is now just the first part of the name from here on.
            out[noutimages]=substr(l_filename,1,strlen(l_filename)-5)
	}
	if (noutimages != nimages) {
	    print("Different number of input images ("//nimages//") and output images ("//noutimages//")")
	    status =1
	    goto clean
	}

        scanfile=""
        delete(tmpfile,ver-, >& "dev$null")
  
    } else {
      print("Error: No output file name(s) were defined.")
      status=1
      goto clean
    }

    if (nbad > 0) {
        print("Error: "//str(nbad)//" of the output file names were bad.")
	status=1
	goto clean
    }

    nbad=0
    i=1
    while (i <= nimages) {
        # check the primary header
	header=in[i]//"[0]"
	imgets(header,"INSTRUMENT")
	instrument=imgets.value
 	if (instrument == "michelle") {
 	  imgets(header,"MPREPARE", >& "dev$null")
 	  if (imgets.value == "0") {
            print("Error: Image "//in[i]//" not MPREPAREd.",
	    status=1
 	    goto clean
 	  }
 	}

	# find the observation mode
	#
        if (instrument == "michelle") {
	    imgets(header,"MODE", >& "dev$null")
	    if (imgets.value == "0") {
	        glogprint( l_logfile, "mipolstack", "status", type="error",
		    errno=131, str="Could not find the MODE from the primary \
		    header.",verbose+)
        	status=status+1
        	goto nextimage
	    }
	    # Change these according to the Michelle "MODE" keywords.
	    # I am not sure whether there are other "non-destructive" (nd) 
	    # modes, or whether these keywords are all correct.
	    modeflag=0
	    if (imgets.value == "chop-nod") modeflag=1
	    if (imgets.value == "ndchop") modeflag=1
	    if (imgets.value == "chop") modeflag=2
	    if (imgets.value == "nod") modeflag=3
	    if (imgets.value == "ndstare") modeflag=4
	    if (imgets.value == "stare") modeflag=4
	    if (modeflag == 0) {
	        glogprint( l_logfile, "mipolstack", "status", type="error",
		    errno=132, str="Unrecognized MODE ("//imgets.value//") \
	            in the primary header.",verbose+)
                status=status+1
                goto nextimage
	    }
        }
        else {
	  print("Error: The instrument used is not MICHELLE")
           status=status+1
           goto nextimage
        }

        # Count the number of extensions

	l_extensions=1
	while (imaccess(in[i]//"["//l_extensions//"]")) {
	    imgets(in[i]//"["//l_extensions//"]","i_naxis")
	    if (modeflag == 1 || modeflag == 2) {
        	if (imgets.value != "3" && l_extensions > 0) {
        	  status=status+1
        	  goto nextimage
        	}
	    }
	    if (modeflag == 3 || modeflag == 4) {
        	if (imgets.value != "2" && l_extensions > 0) {
        	    status=status+1
        	    goto nextimage
        	}
	    }
	    l_extensions=l_extensions+1
	}

	j=l_extensions-1
	print("Number of extensions is "//j)

        if (j < 1) {
	    print ("No data extensions in file "//in[i]//".")
            status = status + 1
            goto nextimage
        }

        # WARNING: modeflag can change within the if (modeflat==?) blocks
        #        therefore, the if-else-if structure SHOULD NOT be used.
        
        if (modeflag == 1) {
            if (8*(j/8) != j) {
              print("Number of extensions for input file "//in[i]//" does not correspond to complete cycles of the polarimetry mode.  Skipping file.")
              goto nextimage
            }
            nnods=j/8
	    if (2*(nnods/2) != nnods) {
	        if (nnods != 1) {
		    print("Number of extensions for input file "//in[i]//" does not correspond to complete nodsets.  Removing last unmatched nod position.")
        	    l_extensions=l_extensions-8
                    if (l_extensions == 0) {
		        print ("No useable data extensions in file "//in[i]//".")
                        goto nextimage
                  }
	        } else
                  modeflag = 2
	    }

# in this initial routine I will NOT check for "BADNOD" flags....everything get stacked up
            if (nnods != 1) {
              for (m=1; m <= nnods; m=m+1) {
                for (j=8*(m-1)+1; j <= 8*m ; j=j+1) {
                  jwp=0
                  k=j-8*(j/8)
                  if (k == 0) k=8
                  if (k == 1 || k == 4) jwp=1
                  if (k == 5 || k == 8) jwp=2
                  if (k == 2 || k == 3) jwp=3
                  if (k == 6 || k == 7) jwp=4
# These are flags for positions 0, 45, 22.5, and 67.5 respectively
                  if (jwp > 0 && jwp < 5) {
                    nwp[jwp]=nwp[jwp]+1
                    filename=out[i]//"_wp"//jwp//"im"//nwp[jwp]
                    imcopy(in[i]//"["//j//"][*,*,3]",filename,ver-, >& "dev$null")
                    if (jwp == 1) {
                       hedit(filename,"WPANGLE",0.0, add+,addonly+,verify-,show-,update+)
                       hedit(filename,"POlANGLE",0.0, add+,addonly+,verify-,show-,update+)
                    }
                    if (jwp == 2) {
                       hedit(filename,"WPANGLE",22.5, add+,addonly+,verify-,show-,update+)
                       hedit(filename,"POlANGLE",45.0, add+,addonly+,verify-,show-,update+)
                    }
                    if (jwp == 3) {
                       hedit(filename,"WPANGLE",45.0, add+,addonly+,verify-,show-,update+)
                       hedit(filename,"POlANGLE",90.0, add+,addonly+,verify-,show-,update+)
                    }
                    if (jwp == 4) {
                       hedit(filename,"WPANGLE",67.5, add+,addonly+,verify-,show-,update+)
                       hedit(filename,"POlANGLE",135.0, add+,addonly+,verify-,show-,update+)
                    }
                  }
                }
              }
            }
	    
        }
        # jump to here if there is a problem

nextimage:
        i=i+1
    }

clean:
    scanfile=""
    delete(tmpfile,ver-, >& "dev$null")
    delete("tmpin*",ver-, >& "dev$null")

exitnow:
    ;

end