procedure nodstack(inimage)
char inimage       {"",prompt="Input T-ReCS or Michelle image"}
char rawpath       {"",prompt="Path for raw images"}
char outprefa      {"a",prompt="Prefix for nod A output image"}
char outprefb      {"b",prompt="Prefix for nod B output image"}
int status         {0,prompt="Exit status: (0=good, >0=bad)"}
struct* scanfile   {"",prompt="Internal use only"}

begin
char l_inputimage,l_rawpath,l_outprefa,l_outprefb
char tmpa,tmpb,instrument,filename,header,nodpos
int modeflag,i,suma,sumb,j

l_inputimage=inimage
l_outprefa=outprefa
l_outprefb=outprefb
l_rawpath=rawpath

# 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=""

if (!imaccess(l_inputimage) && !imaccess(l_rawpath//l_inputimage)) {
  print("Error -- input file "//l_inputimage//" was not found.")
  status=1
  bye
}

if (imaccess(l_outprefa//l_inputimage) || imaccess(l_outprefb//l_inputimage)) {
  print("Error -- output file "//l_outprefa//l_inputimage//" or ")
  print("        "//l_outpreb//l_inputimage//" already exists.")
  status=1
  bye
}

if (l_rawpath=="" || l_rawpath==" ") 
    filename=l_inputimage
else
    filename=l_rawpath//l_inputimage    

tmpa=mktemp("tmpa")
tmpb=mktemp("tmpb")

header=filename//"[0]"
imgets(header,"INSTRUME")
instrument=imgets.value
print("instruemnt = "//instrument)
if (instrument != "michelle" && instrument != "TReCS") {
    print("Error -- instrument keyword "//instrument//" is not recognised.")
    status=1
    goto clean
}

if (instrument == "michelle") {
    imgets(header,"MPREPARE", >& "dev$null")
    if (imgets.value == "0") {
        print("Warning -- file "//l_inputimage//" has not been MPREPAREd.")
        mprepare(l_inputimage,rawpath=l_rawpath,outprefix="m",ver-, >& "dev$null")
        if (mprepare.status != 0) {
            print("Error -- MPREPARE failed for file "//l_inputimage//".")
            status=1
            goto clean
        }
        l_rawpath=""
        filename="m"//l_inputimage
    }
}

if (instrument == "TReCS") {
    imgets(header,"TPREPARE", >& "dev$null")
    if (imgets.value == "0") {
        print("Warning -- file "//l_inputimage//" has not been TPREPAREd.")
        tprepare(l_inputimage,rawpath=l_rawpath,outimage="t"//l_inputimage,ver-, >& "dev$null")
        if (tprepare.status != 0) {
            print("Error -- TPREPARE failed for file "//l_inputimage//".")
            status=1
            goto clean
        }
        l_rawpath=""
        filename="t"//l_inputimage
    }
}

if (instrument == "TReCS") {
    imgets(header,"OBSMODE", >& "dev$null")
    if (imgets.value == "0") {
        print("Error -- obsmode could not be determined.")
        status=1
        goto clean
    }
    modeflag=0
    if (imgets.value == "chop-nod") modeflag=1
    if (imgets.value == "chop") modeflag=2
    if (imgets.value == "nod") modeflag=3
    if (imgets.value == "stare") modeflag=4
    if (modeflag == 0) {
        print("Error -- unknown OBSMODE value.")
        status=status+1
        goto clean
    }
    if (modeflag == 4 || modeflag == 2) {
        print("Error -- the observation is chop or stare mode.")
        status=status+1
        goto clean
    }
}
else {
    imgets(header,"MODE", >& "dev$null")
    if (imgets.value == "0") {
        print("Error -- mode could not be determined.")
        status=1
        goto clean
    }
    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) {
        print("Error -- unknown OBSMODE value.")
        status=status+1
        goto clean
    }
    if (modeflag == 4 || modeflag == 2) {
        print("Error -- the observation is chop or stare mode.")
        status=status+1
        goto clean
    }
}

suma=0
sumb=0

for (i=1; i < 1000; i=i+1) {
  if (imaccess(filename//"["//i//"]")) {
        if (instrument == "TReCS") {
            header=filename//"["//str(i)//"]"
            imgets(header,"nod")
            nodpos=imgets.value
            header=filename//"["//str(i)//"][*,*,3]"
            if (nodpos == "A") {
                if (suma == 0) {
                    imcopy(header,tmpa,ver-, >& "dev$null")
                }
                else {
                    imarith(tmpa,"+",header,tmpa,ver-, >& "dev$null")
                }
                suma=suma+1
            }
            if (nodpos == "B") {
                if (sumb == 0) {
                    imcopy(header,tmpb,ver-, >& "dev$null")
                }
                else {
                    imarith(tmpb,"+",header,tmpb,ver-, >& "dev$null")
                }
                sumb=sumb+1
            }
        }
        else {
            if (modeflag == 1) 
                header=filename//"["//str(i)//"][*,*,3]"
            else
                header=filename//"["//str(i)//"][*,*]"
            j=i-4*(i/4)
            if (j == 0) j=4
            if (j == 1 || j == 4) {
                if (suma == 0) {
                    imcopy(header,tmpa,ver-, >& "dev$null")
                }
                else {
                    imarith(tmpa,"+",header,tmpa,ver-, >& "dev$null")
                }
                suma=suma+1
            }
            else {
                if (sumb == 0) {
                    imcopy(header,tmpb,ver-, >& "dev$null")
                }
                else {
                    imarith(tmpb,"+",header,tmpb,ver-, >& "dev$null")
                }
                sumb=sumb+1
            }
        }
  }
  else {
    if (i == 1) {
        print("Error -- no data extensions in file "//l_inputimage)
        status=1
        goto clean
    }
    goto finish
  }
}

finish:

  header=filename//"[0]"

  if (suma == 0) {
    print("Error -- no NOD A frames found.")
    status=1
  }
  else {
    imarith(tmpa,"/",suma,tmpa)
    wmef(tmpa,l_outprefa//l_inputimage,extname="SCI",phu=header,verb-, >& "dev$null")
  }
  if (sumb == 0) {
    print("Error -- no NOD B frames found.")
    status=1
  }
  else {
    imarith(tmpb,"/",sumb,tmpb)
    wmef(tmpb,l_outprefb//l_inputimage,extname="SCI",phu=header,verb-, >& "dev$null")
  }

clean:

  imdelete(tmpa,ver-, >& "dev$null")
  imdelete(tmpb,ver-, >& "dev$null")

end