misc_version=1.05

# get_image					- Added by C. Edwards 6/2009
# get_map -- added by C.Edwards 9/3/09
# available_maps -- added by C.Edwards 9/3/09
# getopt        - taken from /u/gorelick/getopt.dv
# opt_usage     - taken from /u/gorelick/getopt.dv
# strsplit			- split a string into a text array based on a delimiter
# uniq          - taken from /u/cedwards/christopher.dvrc, formerly called guniq()
# weight_array	- taken from /u/phil/prc_dvrc
# k_phil	- taken from /u/knowicki/dvrc/k_phil.dvrc
# k_fill	- taken from /u/phil/k_fill.dvrc
# thmproc_find  - taken from /u/knowicki/dvrc/tools.dvrc
# mars_bin      - taken from /u/bandfield/josh.dvrc
# map_trans     - taken from /u/bandfiled/josh.dvrc
# geo_trans     - taken from /themis/lib/themis_dvrcjlb
# p_surf        - taken from /themis/lib/themis_dvrcjlb
# read_vm       - taken from /u/bandfield/josh.dvrc, formerly called readvm()
# write_vm      - taken from /u/bandfield/josh.dvrc, formerly called writevm2()
# make_mtesx	- taken from /mtes/lib/mt_dvrc_cal
# make_tesx	- taken from /tes/lib/tes_dvrc_cal
# make_lab1x	- taken from /u/phil/prc_dvrc
# make_themisx	- taken from /u/phil/prc_dvrc
# lsubs         - Added by J. Bandfield 1/2008
# deconstruct   - Added by J. Bandfield 1/2008
# rm_quote      - Added by PRC 11/2010
# bs						- created by C.Edwards 12/10/2010

#1.04 added rm_quote 

define get_image(instrument,type) {

#modified 10/7/09 to query for crism data.
#modified 5/5/10 to take alternate data types

	if($ARGC==0) {
		printf("\nFind the desired data based on image id/instrument/type\n")
		printf("Returns an url formatted for the defined query\n\n")
		printf("$1 = image id\n\n")
		printf("Optional:\n")
		printf("instrument = \t\t\"ctx\",\"hirise\",\"hrsc\",\"moc\"\n\t\t\t\"viking\",\"anaglyph\",\"themis\",\"crism\" (Default = \"themis\")\n\n")
		printf("type (if hirise) = \t\"RGB\",\"RED\" (Default = \"RED\")\n")
		printf("type (if themis) = \t\"BWS\",\"DCS\",\"D875\",\"D964\",\"D642\"\n\t\t\t\"PBT\",\"BTR\",\"RDR\",\"ABR\" (Default = \"RDR\")\n")
		printf("type (if crism) = \t\"TRDR-LBL\",\"TRDR-IMG\",\"DDR-LBL\",\"DDR-IMG\",\"MAF\"\n")
		printf("\t\t\t\"HYD\",\"PHY\",\"ICE\",\"TRU\",\"IRA\",\"FEM\" (Default = \"TRDR-LBL\")\n")
		printf("\t\tNOTE:\tThe following crism types require a \"_RGB\", \"_HSV\", or\" _BW\" as a suffix\n")
		printf("\t\t\t\"AL-OH\",\"FE-OH\",\"FERRIC\",\"H2O-1\",\"H2O-2\",\"H2O-3\",\"HCP\"\n")
		printf("\t\t\t\"HYD-SI\",\"LCP\",\"OLIVINE\",\"PHYLLO\",\"SULFATE\",\"ZEOLITE\"\n\n")
		printf("type (if hrsc) = \t\"DT4\", \"ND3\", \"ND4\", \"SR3\" (Default = \"ND3\")\n\n")
		printf("\t\tNOTE:\ttype is only used for THEMIS/CRISM/hrsc/hirise data\n\n")
		printf("c.edwards 3/09\n")
		return(null)
	}
	
	id=$1
	if(HasValue(instrument)==0) instrument="themis"
	if(HasValue(type)==0 && instrument=="themis") type="RDR"
	if(HasValue(type)==0 && instrument=="crism") type="TRDR-LBL"
	if(HasValue(type) && instrument=="hrsc") {
		id=id+"_"+type
		type=""
	}
	if(HasValue(type) &&  instrument=="hirise") { 
		id=id+"_"+type
		type=""
	}
	if(HasValue(type)==0) type=""

  if(HasValue(devserver)==0) { 
		server="http://app-davinci.mars.asu.edu/StampServer-1.2/ImageServer?"
	}

	ver=version()[18:]
	ver="davinci "+ver

	url=sprintf("%sinstrument=%s&id=%s&imageType=%s&version=%s",server,instrument,id,type,ver)
	return(url)
}


define available_maps(format,search,mapname) {

	#populate default options
	if(HasValue(format)==0) format=1
	if(HasValue(search)==0) search=""
	if(HasValue(mapname)==0) mapname=""
			
	#get the mapserver xml and ghetto parse it
	verbose=0
	xml=read_lines("http://ms.mars.asu.edu/?SERVICE=WMS&REQUEST=GetCapabilities")
	verbose=3
	names=grep(xml,"<Name>")[,2:]
	titles=grep(xml,"<Title>")[,3:]
	lines=clone(0,1,length(names),1)

	namestart=strstr(names,"<Name>")+6
	nameend=strstr(names,"</Name>")-1
	titlestart=strstr(titles,"<Title>")+7
	titleend=strstr(titles,"</Title>")-1
	
	if(mapname!="")	oldnames=names

	#loop through the names to find the beginning and ends
	for(i=1;i<=length(names);i+=1) {
		names[,i]=names[namestart[,i]:nameend[,i],i]			
		titles[,i]=titles[titlestart[,i]:titleend[,i],i]
	}

	if(mapname!="") {	
		for(i=1;i<=length(names);i+=1) {
			if(names[,i]==mapname) { 
				xmlmapname=oldnames[,i]
				xmlnextmapname=oldnames[,i+1]
			}
		}
		for(i=1;i<=length(xml);i+=1) {
			if(xmlmapname==xml[,i]) mapline=i
			if(xmlnextmapname==xml[,i]) nextmapline=i	
		}
		numeric=grep(xml[,mapline:nextmapline],"<Keyword>numeric")
		if(HasValue(numeric)) {
			dl_type="vicar"
		} else {
			dl_type="png"
		}
		return(dl_type)
	}
		
	for(i=1;i<=length(names);i+=1) {
		for(j=1;j<=length(xml);j+=1) {
			if(names[,i]==xml[,j]) {
				lines[,i]=j
			}
		}
	}

	#search option
	titles2=""
	if(search!="") {
		names2=grep(names,search)
		for(i=1;i<=length(names2);i+=1) {
			for(j=1;j<=length(names);j+=1) {
				if(names2[,i]==names[,j]) {
					titles2=cat(titles2,titles[,j],axis=y)
				}
			}
		}
		titles=titles2[,2:]
		names=names2
	}

	#sort alphabetically
	titles=sort(titles,by=names)
	names=sort(names)	

	#First output option - table view
	if(format==1) {
		names=cat("Map Layer Name",names,axis=y)
		titles=cat("Map Title",titles,axis=y)

		titlemaxlen=int(max(strlen(titles)))
		namemaxlen=int(max(strlen(names)))
  	for (i = 1 ; i <= length(names) ; i++) {
			spaces=" "
			titlecurlen=strlen(titles[,i])
			namecurlen=strlen(names[,i])

			for(j=titlecurlen;j<=titlemaxlen;j+=1) {
				titles[,i]=cat(titles[,i]," ",axis=x)
			}
			for(j=namecurlen;j<=namemaxlen-1;j+=1) {
				spaces=cat(spaces," ",axis=x)
			}
		
			if(i==1) {
				breaks="-"
				for(j=1;j<=titlemaxlen+namemaxlen+9;j+=1) {
					breaks=cat(breaks,"-",axis=x)
				}
				printf("%s\n",breaks)
				printf('| %s | %s %s |\n', ""+titles[,i], ""+names[,i],spaces);
				printf("%s\n",breaks)
			} else if (i==length(names)) {
				printf('| %s | "%s"%s|\n', ""+titles[,i], ""+names[,i],spaces);
				printf("%s\n",breaks)
			} else {
				printf('| %s | "%s"%s|\n', ""+titles[,i], ""+names[,i],spaces);
			}
		}
		printf("\nTry format=2 for a different view, or format=3 to return a structure\n")
		printf("Also search by search=\"THEMIS\"\n")
	}

	#second output option - list view
	if(format==2) {
	  str = "  ";
 		for (i = 1 ; i <= length(names) ; i++) {
			name=names[,i]
    	s = sprintf('"%s", ', ""+name);
    	if (length(str+s) > 100) {
      	printf("%s\n", str)
      	str = "    " + s
    	} else {
      	str = str + "  " + s
    	}
  	}
  	printf("%s\n", str)	
	}
	
	#third output option (return the values as a struct)
	if(format==3) {
		out={}
		out.names=names
		out.titles=titles
		return(out)
	}
}
	


define get_map(lat,lon,clat,clon,lpo,spo,ppd,proj,map,trim,size,struct,v){

	#help
	if((HasValue(lat)==0 || HasValue(lon)==0) && $ARGC==0) {
	  printf("\nGet Map data from ms.mars.asu.edu\n")
 		printf("Will download map data and load it as a variable\n\n")
		printf("Can either take an ISIS cube or supplied lat/lon range\n")
		printf("\t$1=ISIS structure (read with load_pds())\n")
		printf("or\n")
		printf("\tlat=latitude range (\"min:max\")\n")
		printf("\tlon=longitude range (\"min:max\")\n")
		printf("\tppd=pixels per degree (default=8)\n\n")
		printf("Options:\n")
		printf("\tmap=ms.mars.asu.edu map name (default=\"MOLA_128ppd_shade_ne\")\n")
		printf("\tproj=projection type:\n")
		printf("\t\t\"SIMPLE_CYLINDRICAL\"\n")
		printf("\tclat=center latitude (Default=average latitude)\n")
		printf("\tclon=center longitude (Default=average longitude)\n")
		printf("\tlpo=line projection offset (Default=0)\n")
		printf("\tspo=sample projection offset (Default=0)\n")
		printf("\ttrim=trim the output data to the ISIS input\n")
		printf("\tsize=pixels per tile (default=256)\n")
		printf("\tv=verbosity (0,1,2)\n")
		printf("\tstruct=return a structure of parameters (0,1)\n")
		printf("\nFor Available Maps See:\n")
		printf("available_maps()\n")
		printf("\ncedwards 3/2009\n")
		return(null)
	}

	#Check to see if we are using an ISIS header
	#make sure we don't define things many time
	if($ARGC>=1) {
		if(type($1)=="STRUCT") {
			 pds=$1
		} else {
			printf("ISIS structure not supplied\n")
			return(null)
		}
	}

	if(HasValue(ppd) && $ARGC>=1) {
		printf("ppd cannot be specified with an ISIS structure\n")
		return(null)
	}
	if(HasValue(proj) && $ARGC>=1) {
		printf("proj cannot be specified with an ISIS structure\n")
		return(null)
	}
	if(HasValue(clat) && $ARGC>=1) {
		printf("clat cannot be specified with an ISIS structure\n")
		return(null)
	}
	if(HasValue(clon) && $ARGC>=1) {
		printf("clon cannot be specified with an ISIS structure\n")
		return(null)
	}
	if(HasValue(lpo) && $ARGC>=1) {
		printf("lpo cannot be specified with an ISIS structure\n")
		return(null)
	}
	if(HasValue(spo) && $ARGC>=1) {
		printf("spo cannot be specified with an ISIS structure\n")
		return(null)
	}
	if(HasValue(lat) && $ARGC>=1) {
		printf("lat cannot be specified with an ISIS structure\n")
		return(null)
	}
	if(HasValue(lon) && $ARGC>=1) {
		printf("clon cannot be specified with an ISIS structure\n")
		return(null)
	}
	if(HasValue(trim) && $ARGC==0) {
		printf("trim cannot be specified without an ISIS structure\n")
		return(null)
	}

	#populate default options 
 	if(HasValue(map)==0) map="MOLA_128ppd_shade_ne"
  if(HasValue(ppd)==0) ppd=8
	if(HasValue(size)==0) size=256
	if(HasValue(v)==0) v=0
	if(HasValue(trim)==0) trim=0
	if(HasValue(proj)==0) proj="SIMPLE_CYLINDRICAL"	
	if(HasValue(struct)==0) struct=1

	#fill out some default values in param structure
	param={}
	param.map=map	
	param.ppd=ppd
	param.size=size
	param.v=v
	param.trim=trim

	#make sure the map you chose is a possible option
	maps=available_maps(format=3).names
	if(HasValue(grep(maps,map+"$"))==0) {
		printf("Map: %s is not a valid map layer\n",map)
		printf("Run available_maps() to see a list of possible maps\n")
		printf("Or search for maps using available_maps(search=\"GRS\")\n")
		return(null)
	}
	
	verbose=0
	#if we are using a user supplied range
	if($ARGC==0) {
		#make sure lat/lon ranges are correct
		if(strstr(lat,":")==0) {
 		 	printf("No latitude range specified (\"min:max\")\n");
 		 	return(null);
		}
		if(strstr(lon,":")==0) {
 		 	printf("No longitude range specified (\"min:max\")\n");
 			return(null);
		}

		#build the latitude range
		latpos=strstr(lat,":")
		param.minlat=atof(lat[:latpos-1])
		param.maxlat=atof(lat[latpos+1:])
		if(param.minlat>param.maxlat) {
 		 	lattmp=param.maxlat
 	 		param.maxlat=param.minlat
 	 		param.minlat=lattmp
		}

		#build longitude range
		lonpos=strstr(lon,":")
		param.minlon=atof(lon[:lonpos-1])
		param.maxlon=atof(lon[lonpos+1:])
		if(param.minlon>param.maxlon) {
	 	 	lontmp=param.maxlon
 		 	param.maxlon=param.minlon
 		 	param.minlon=lontmp
		}
	
		#see if the user has defined clat/clon/lpo/spo and fill out defaults if not
		if(HasValue(clat)==0) clat=(param.minlat+param.maxlat)/2.0
		if(HasValue(clon)==0) clon=(param.minlon+param.maxlon)/2.0

		#these may not be right....
		if(HasValue(lpo)==0) lpo=(-1)*param.maxlat*ppd
		if(HasValue(spo)==0) spo=param.minlon*ppd
#		if(HasValue(spo)==0) spo=(param.minlon-clon)*ppd*cos(param.maxlat)

	} else {

		#ISIS case
		#begin extracting values form the image_map_projection keyword
		struct=pds.qube.image_map_projection
		param.minlat=struct.minimum_latitude
		param.maxlat=struct.maximum_latitude
		param.minlon=struct.westernmost_longitude
		param.maxlon=struct.easternmost_longitude

		#make sure things go from small to large
		if(param.maxlon<param.minlon) {
			param.minlon=param.maxlon
			param.maxlon=struct.westernmost_longitude
		}
	
		#extract more projection parameters
		param.ppd=struct.map_resolution
		spo=struct.sample_projection_offset
		lpo=struct.line_projection_offset

		#make values if structure does not contain them
		if(HasValue(struct.center_latitude)==0) {
			clat=(param.maxlat-param.minlat)/2.0
		} else {
			clat=struct.center_latitude
		}
		if(HasValue(struct.center_longitude)==0) {
			clat=(param.maxlon-param.minlon)/2.0
		} else {
		clon=struct.center_longitude
		}

		#extract the projection (without the second set of quotes)
		proj=struct.map_projection_type[2:length(struct.map_projection_type)[1]-1]
	}


	#change the projection from ISIS notation to JMARS notation
	param.projection=proj
	if(proj=="SIMPLE_CYLINDRICAL") {
		param.jproj="JMARS:1"
#	} else if (proj=="SINUSOIDAL") {
#		param.jproj="JMARS:0"
#	} else if (proj="POLAR_STEREOGRAPHIC") {
#		param.jproj="JMARS:2"
	} else {
		printf("Error: Incorrect Projection Supplied\n")
		return(null)
	}

	verbose=0

	#fill out other paramters from ISIS or user defined ranges
	param.lpo=lpo
	param.spo=spo
	param.clat=clat
	param.clon=clon

	#param.dl_type=available_maps(mapname=map)
	verbose=0
	param.dl_type="vicar"

	#calculate the size from lat/lon/ppd and test chunk dim
	param.x=int(ceil(param.ppd*abs(param.minlon-param.maxlon)))
	param.y=int(ceil(param.ppd*abs(param.minlat-param.maxlat)))
	url="http://ms.mars.asu.edu/?SERVICE=WMS&REQUEST=GetMap&FORMAT=image/"+param.dl_type+"&WIDTH=128&HEIGHT=128&SRS=JMARS:1&BBOX=0,0,5,5&STYLES=&VERSION=1.1.1&LAYERS="+map
	chunk=read(url)
	param.z=dim(chunk)[3]
	if(param.z==4) param.z=3
	if(param.z==2) param.z=1
	if(v==2) printf("Size: %ix%ix%i\n",param.x,param.y,param.z)


	#create the output data of type and org with the test chunk	
	data=org(format(clone(0,param.x,param.y,param.z),format(chunk)),org(chunk))
	count=0

	#start building up the chunks
	for(i=1;i<=param.x;i+=param.size) {
	
		#define the lon range and the size of the box
		lon1=(i-1)/param.ppd+param.minlon
		lon2=(i+param.size)/param.ppd+param.minlon
		if(lon2>param.maxlon) lon2=param.maxlon
		xsize=int(param.ppd*(lon2-lon1))
		lpoi=(i-1)+param.lpo

		for(j=1;j<=param.y;j+=param.size) {

			#define the lat range and the size of the box
			lat1=param.maxlat-(j-1)/param.ppd
			lat2=param.maxlat-(j+param.size)/param.ppd
			if(lat2<param.minlat) lat2=param.minlat
			ysize=int(param.ppd*(lat1-lat2))
			spoi=(j-1)+param.spo

			#verbose options
			if(v>=1) {
				count+=1
				printf("%i of %i\n",count,ceil(param.x/float(param.size))*ceil(param.y/float(param.size)))
			}
			if(v==2) printf("lon=\"%f:%f\", lat=\"%f:%f\", size=\"%i,%i\", range=\"%i:%i,%i:%i\", offset=\"%i,%i\"\n",lon1,lon2,lat1,lat2,xsize,ysize,i,i+xsize-1,j,j+ysize-1,lpoi,spoi)
	
			#build the mapserver url
			if(param.jproj=="JMARS:1") {

				#SIMPLE_CYLINDRICAL CASE
				url="http://ms.mars.asu.edu/?SERVICE=WMS&REQUEST=GetMap&FORMAT=image/"+param.dl_type+"&WIDTH="+xsize+"&HEIGHT="+ysize+"&SRS="+param.jproj+"&BBOX="+lon1+","+lat2+","+lon2+","+lat1+"&STYLES=&VERSION=1.1.1&LAYERS="+map

#			} else if (param.jproj=="JMARS:0" || param.jproj=="JMARS:2") {
#				#POLAR_STEROGRAPHIC AND SINUSOIDAL CASES
#				url="http://ms.mars.asu.edu/?SERVICE=WMS&REQUEST=GetMap&FORMAT=image/"+param.dl_type+"&WIDTH="+xsize+"&HEIGHT="+ysize+"&SRS="+param.jproj+","+param.clon+","+param.clat+","+spoi+","+lpoi+","+param.ppd+"&BBOX="+lon1+","+lat2+","+lon2+","+lat1+"&STYLES=&VERSION=1.1.1&LAYERS="+map

			} else {
				printf("Error: Wrong projection\n")
				return(null)
			}

			#get the data and fill in the chunks
			data[i:i+xsize-1,j:j+ysize-1,:zdim]=read(url)[,,:zdim]
		}
	}
	verbose=3

	if(type(data)=="byte") {			
		param.ignore=0
	} else {
		param.ignore=-32768
	}

	if(trim==1) {
		data[where pds.qube.data==min(pds.qube.data)]=param.ignore
	}

	#if we have a structure then fill it and return
	if(struct==1) {
		param.ul=param.maxlat//param.minlon
		param.lr=param.minlat//param.maxlon
		param.ur=param.maxlat//param.maxlon
		param.ll=param.minlat//param.minlon
		param.data=data
		return(param)
	} else {	
		return(data)
	}
}



define getopt() {

  if($ARGC==0) { 
    printf("\nThis function takes a structure containing the available options")
    printf("in the expected order, and fills in values\n")
    printf("Will also grab the command line options and fill in values (use two dashes)\n\n")
    printf("n.gorelick\n")
    return(null)
  }

  args = $1
  args_copy = args;
  if ($ARGC == 2) {
    opt = $2
  } else {
    opt = {}
  }
  out = opt
  opt_s = {}
  for (i = length(args) ; i > 0 ; i--) {
    arg = args[i]
    if (strstr(arg, "=")) {
      pos = strstr(arg, "=")
      name = arg[1:pos-1]
      value = arg[pos+1:]
      while(name[1] == '-') name = name[2:]
      if (HasValue(get_struct(out, name))) {
        form=format(get_struct(out,name))
	remove_struct(out,name)
      }
      if(HasValue(form)==0) {
 	add_struct(out, name=name, value=value);
      } else {
        if(form=="STRING") {
          add_struct(out, name=name, value=value);
        } else {
          add_struct(out, name=name, value=format(atof(value),format=form));
        }
      }
      add_struct(opt_s, name=name, value=1);
      remove_struct(args_copy, pos=i);
    }
  }
  args = args_copy
  count = 1;
  for (i = 1 ; i < length(opt) ; i++) {
    name = get_struct_key(opt, i)
    if (count <= length(args) && get_struct(opt_s, name) != 1) {
      remove_struct(out, name);
      add_struct(out, name=name, value=args[count]);
      count++;
    }
  }
  return(out)
}

define opt_usage() {

   if($ARGC==0) { 
    printf("\nThis function prints the available options to the screen\n")
    printf("Used in conjunction with getopt()\n")
    printf("n.gorelick\n")
    return(null)
  }

  opt = $1
  str = "  ";
  for (i = 0 ; i < length(opt) ; i++) {
    name = get_struct_key(opt, i+1)
    value = opt[i+1];

    s = sprintf('%s="%s"', name, "" + value);
    if (length(str+s) > 72) {
      printf("%s\n", str)
      str = "    " + s
    } else {
      str = str + "  " + s
    }
  }
  printf("%s\n", str)
}

define strsplit(delim) {
	if($ARGC==0) {
	  printf("\nSplit a string or string array into an array of strings or structure of strings by a delimiter\n\n")
    printf("$1=string\n")
		printf("delim=the splitting delimeter (Default is \",\")\n\n")
    printf("c.edwards 05-06-10\n")
    return(null)
	}
  string=$1
  if(HasValue(delim)==0) {
    delim=","
  }

	if(strstr(string,delim)==0) {
 	   printf("No delimiter was found\n")
 	   return(null)
 	 }

	if(type(string)=="STRING") {
	  out=text()
 	 	while(strstr(string,delim)!=0) {
 	   	pos=strstr(string,delim)
 	   	out=cat(out,string[:pos-1],axis=y)
 	   	string=string[pos+1:]
 	 	}
 	  out=cat(out,string,axis=y)
	} else if (type(string)=="TEXT"){

		struct={}
		for(i=1;i<=length(string);i+=1) {

		  out=text()
 		 	while(strstr(string[,i],delim)!=0) {
 	   		pos=strstr(string[,i],delim)
 	   		out=cat(out,string[:pos-1,i],axis=y)
 	   		string[,i]=string[pos+1:,i]
	 		}

 	 		out=cat(out,string[,i],axis=y)
			insert_struct(struct,name=i+"",value=out)
		}

		out=struct

	} else {
		printf("Invalid data type\n")
		return(null)
	}
  return(out)
}



define uniq() {

  if($ARGC==0) { 
    printf("\nGhetto find the uniq numbers in a one dimensional array\n")
    printf("A sorted one dimensional array of the same input format and direction will be returned\n\n")
    printf("$1=data array\n\n")
    printf("c.edwards 03-22-07\n")
    return(null)
  }
  
  verbose=0
  #get the uniq answers from one dimensional data
  data=$1
  dim=dim(data)
  max=max(dim)
  dim[where dim!=max]=0
  struct={}  
   
  for(dir=1;dir<=3;dir+=1) {
    if(dim[dir]==max) break;
  }

  if(dir==1) data=translate(data,x,z)
  if(dir==2) data=translate(data,y,z)
  
  for(i=1;i<=max;i+=1) {
    add_struct(struct,name=""+data[,,i])
  }
  
  newdata=clone(format(0,format(data)),1,1,length(struct))
  for(i=1;i<=length(struct);i+=1) {
    newdata[,,i]=format(atof(get_struct_key(struct,i)),format(data))
  }
  
  if(dir==1) newdata=translate(newdata,z,x)
  if(dir==2) newdata=translate(newdata,z,y)

  verbose=3
  return(sort(newdata))

}



define weight_array() {

if ($ARGC == 0) {
	printf("make weighting function for use in thm.kfill\n")
	printf(" \n")
	printf("usage:	wt = weight_array(radius, power) \n")
	printf(" radius in pixels\n")
	printf(" power = 0 - linear weighting with distance\n")
	printf(" power > 0 - weight = 1/(distance * power)\n")
	printf(" \n")
	return(0)
}
# p.christensen
radius = $1
power = $2

nxy = int((radius * 2) + 1)
center = radius + 1
R = sqrt(radius^2. + radius^2.)
w = create(nxy, nxy, 1, format=float)

for(i=1; i<= nxy; i+=1) {
	for(j=1; j<= nxy; j+=1) {
		# compute distance - then weight by user-input power
		dist = sqrt( (center - i) ^2. + (center - j) ^2. )

		if(power == 0) {
			# do linear case
			w[i,j] = (radius - dist + 1.) / radius
			if(w[i,j] < 0.) {
				w[i,j] = 0.
			}
	
		} else {
			w[i,j] = 1./(dist^power )
		}
	}
}

# set center point to 1.
w[center,center] = 1.
return(w)
}



define k_phil(width,order,iterations) {

  if ($ARGC == 0) {
    printf("Fills the holes in pictures with iterpolated data\n\n")
    printf("USE: b = k_phil(a)\n\n")
    printf("options width, order, and iterations\n")
    printf("width changes the size of the weighting kernel\n")
    printf("order changes the dependence on r of the weighting kernel\n")
    printf("iterations sets the number of iterations the program should run\n\n")
    printf("example: b = k_phil(a,width=15,order=3,iterations=5)\n")
    return(null)
  }

  pic=$1 #the picture ready to be fabricated for monetary gain

  wth=6
  ord=2
  its=0

  if(HasValue(width)) {
    wth=width
  }

  if(HasValue(order)) {
    ord=order
  }

  if(HasValue(iterations)) {
    its=iterations
  }

  w=weight_array(wth,ord) # our weighting array

  if(dim(pic)[2]>500) {
    l=sum(pic,axis=x)
    p_top=l[,:250]
    l_dim=dim(l)
    p_bot=l[,l_dim[2]-250:]
    p_top=translate(p_top,from=y,to=y,flip=1)
    p_top=250-minchan(p_top)+1
    p_bot=l_dim[2]-250+minchan(p_bot)

    top=pic[,:p_top]
    bot=pic[,p_bot:]
    pic=pic[,p_top+1:p_bot-1]
  }

  pic=byte(thm.kfill(pic,w,its))

  if(dim(pic)[2]>500) {
    pic=cat(top,pic,bot,axis=y)
  }

  return(pic)
}




define k_fill(width,order,iterations) {

  verbose=0
  if ($ARGC == 0) {
    printf("Fills the holes in pictures with iterpolated data\n\n")
    printf("USE: b = k_fill(a)\n\n")
    printf("options width, order, and iterations\n")
    printf("width changes the size of the weighting kernel\n")
    printf("order changes the dependence on r of the weighting kernel\n")
    printf("iterations sets the number of iterations the program should run\n\n")
    printf("example: b = k_fill(a,width=15,order=3,iterations=5)\n")
    return(null)
  }

  pic=$1 #the image

  wth=6
  ord=2
  its=0

  if(HasValue(width)) {
    wth=width
  }

  if(HasValue(order)) {
    ord=order
  }

  if(HasValue(iterations)) {
    its=iterations
  }

  w=weight_array(wth,ord) # our weighting array

  pic=byte(thm.kfill(pic,w,its))

  return(pic)
}

define thmproc_find(id,batch) {

  if ($ARGC == 0 && HasValue(id) == 0 && HasValue(batch) == 0) {
    printf("thmproc_find() - Fri Jan  6 17:20:17 MST 2006\n\n")
    printf("returns all internal urls to jobs run in thmproc\n")
    printf("USE: thmproc_find(id,batch)\n\n")
    printf("     thmproc_find(id=\"I01002002\")\n")
    printf("     thmproc_find(batch=\"1314\")\n\n")
    printf("this is an either/or kind of thing.\n\n")
    printf("Requires the user to have database access and SELECT access on horus\n")
    return(null)
  }

  verbose=0
  if(HasValue(id)) {
		args=syscall("themis_db_access --arglist")[,1]
		list=syscall("mysql "+args+" -D horus -B -s -e 'select internal_url,status from cache_item where product_id = \""+id+"\";'")
		
    list=list[8:]
    l = length(list)
    for(i=1;i<=l;i+=1) {
      printf(list[,i]+"\n")
    }

  } else if(HasValue(batch)) {
		args=syscall("themis_db_access --arglist")[,1]
		list=syscall("mysql "+args+" -D horus -B -s -e 'select ci.internal_url,ci.status from cache_item ci, job_component jc, batch_job bj where ci.cache_id = jc.cache_id and jc.job_id = bj.job_id and bj.batch_id = \""+batch+"\";'")

    list=list[8:]
    l = length(list)
    for(i=1;i<=l;i+=1) {
      printf(list[,i]+"\n")
    }

  }
  verbose=3
  return(list)
}



define mars_bin(nofile,null,polar) {

if ($ARGC == 0) {
	printf (" \n")
	printf (" bin tes data \n")
	printf (" $1 = upper lon \n")
	printf (" $2 = lower lon \n")
	printf (" $3 = upper lat \n")
	printf (" $4 = lower lat \n")
	printf (" $5 = resolution (pixels per degree) \n")
	printf (" $6 = vanilla filename or vanilla data directory \n")
	printf (" $7 = binned fields string if directory is specified (optional) \n")
	printf (" $8 = select string if directory is specified (optional) \n")
	printf (" nofile = 1; Use array instead of file for source data.\n")
	return(0)

}


if ($ARGC == 8) {

system("vanilla -select '" + $8 + "' -fields 'longitude latitude " + $7 + "' " + $6 + " | tail +2 > temp1.txt")

van=ascii(filename="temp1.txt", format=double)

}

if (HasValue(nofile)==0) {

# If user supplies vanilla emissivity file

if ($ARGC == 6) {
van=ascii(filename=$6, format=float)

}
}

if (HasValue(nofile)==1 && $ARGC == 6) {

	van=$6

}

if ($ARGC == 3) {

        resolution=$2
        radius=(90-abs($1))*resolution
        van=$3
}


# Get dimensions of the vanilla file

vecs=dim(van)
nlines=vecs[2,,]
nsamples=vecs[1,,]-2

if (HasValue(polar)==0) {

        polar=0
}

if (polar == 0) {
        # Get dimensions of lat and lon box

        ulx=$1
        lrx=$2
        uly=$3
        lry=$4
        resolution=$5

        # translate lat lon samples to data cube range

        lon=int((van[1,,]-lrx)*resolution)+1
        lat=int((van[2,,]-lry)*resolution)+1

        # Figure out the size of the box

        lon_range=int((ulx-lrx)*resolution)
        lat_range=int((uly-lry)*resolution)
}

if (polar==1) {

        # Get dimensions of the lat and lon box

        lon_range=int(radius*2)
        lat_range=int(radius*2)

        # translate lat lon samples to data cube range (lat is y lon is x)

        lon=int(radius+(sind(van[1])*(90-van[2])*resolution))+1
        lat=int(radius-(cosd(van[1])*(90-van[2])*resolution))+1
}

lon [ where (lon == 0.) ] = -1
lat [ where (lat == 0.) ] = -1

# create number of obs array and data array

num_obs=create(x=lon_range, y=lat_range, z=1, format=float, start=0, step=0)

data_cube=create(x=lon_range, y=lat_range, z=nsamples, format=double, start=0, step=0)


# Start looping through the vanilla data and adding to num_obs and data_cube arrays


if (nsamples > 1) {

	for (i=1; i<=nlines; i+=1) {

		num_obs[lon[1,i,1],lat[1,i,1],1] = num_obs[lon[1,i,1],lat[1,i,1],1] + 1

		data_cube[lon[1,i,1],lat[1,i,1],] = data_cube[lon[1,i,1],lat[1,i,1],] + translate(van[3:(nsamples+2),i,1],from=x, to=z)

	}
}

if (nsamples == 1) {

	for (i=1; i<=nlines; i+=1) {

		num_obs[lon[1,i,1],lat[1,i,1],1] = num_obs[lon[1,i,1],lat[1,i,1],1] + 1

		data_cube[lon[1,i,1],lat[1,i,1],] = data_cube[lon[1,i,1],lat[1,i,1],] + van[3,i,1]

	}
}

# divide data_cube by number of observations for each sample

data_cube=data_cube/num_obs

if (polar ==0) {
        # Flip each image to west longitude

        data_cube=translate(data_cube, from=x, to=x, flip=1)
        num_obs=translate(num_obs, from=x, to=x, flip=1)

        data_cube=translate(data_cube, from=y, to=y, flip=1)
        num_obs=translate(num_obs, from=y, to=y, flip=1)
}
# create data output structure

binout=struct(cube,nobs)

binout.cube=data_cube
binout.nobs=num_obs

return(binout)



# clean up if necessary

if ($ARGC == 8) {
system("rm temp1.txt")
}

}



define map_trans() {

if ($ARGC == 0) {
	printf (" \n")
	printf (" Translate image coordinates to map coordinates \n")
	printf (" First pixel is x=0 y=0 \n")
	printf (" $1 = image x \n")
	printf (" $2 = image y \n")
	printf (" $3 = image resolution (ppd) (optional, default=1)\n")
	printf (" $4 = ulx (optional, default=180) \n")
	printf (" $5 = uly (optional, default=90) \n")
	return(0)

}


if ($ARGC == 2) {

	ulx=float(180)
	uly=float(90)
	image_res=float(1)

}

if ($ARGC == 3) {

	ulx=float(180)
	uly=float(90)
	image_res=float($3)

}


if ($ARGC == 5) {

	ulx=float($4)
	uly=float($5)
	image_res=float($3)

}

image_x=float($1)
image_y=float($2)

lower_x=ulx-((image_x+1)/image_res)
upper_x=ulx-((image_x)/image_res)

lower_y=uly-((image_y+1)/image_res)
upper_y=uly-((image_y)/image_res)

if (upper_x <= 0) {

	upper_x=360+upper_x
	lower_x=360+lower_x

}

printf ("Latitude is " + lower_y + " to " + upper_y + " \n")
printf ("Longitude is " + lower_x + " to " + upper_x + " \n")
return(0)

}



define geo_trans() {

if ($ARGC == 0) {
        printf (" \n")
        printf (" Translate map coordinates to image coordinates \n")
        printf (" First pixel is x=1 y=1 \n")
        printf (" $1 = latitude \n")
        printf (" $2 = longitude \n")
        printf (" $3 = image resolution (ppd)\n")
        printf (" $4 = center longitude of image (optional, default is 0)\n")
        return(0)

}


lat=$1
lon=$2
image_res=float($3)

if ($ARGC == 3) {

        center=0
}

if ($ARGC == 4) {

        center=$4
}

image_y=int((90-lat)*image_res+1)
image_x=int((center+180-lon)*image_res)+1

if ((center+180-lon) <= 0) {

        image_x=int(360*image_res)+image_x-1
}

return(image_x//image_y)

}




define p_surf() {

if ($ARGC == 0) {

	printf (" \n")
	printf (" Get surface pressure from elevation and sclk time\n")
	printf (" $1 = elevation (km) \n")
	printf (" $2 = sclk time \n")
	return(0)

}

# This is borrowed by code originally written by Hugh Kieffer in f77, modified by
# J. Pearl, and the modified by M. Kaelbarer and adapted to f90.

# J. Pearl used a scale factor of 0.8979 for TES surface pressure to stabilize
# PT retrievals.  This is not done here.

# Get Julian Date

date=4239.5 + (double($2)/86400)

# Set length of Martian year

marsyear=686.98

# Set Julian date for start of curve

start_julian=3395.49

# Find fraction of year

frac=((date-start_julian)%marsyear)/marsyear

# Set amplitude and phase of harmonic terms for normalized relative pressure from
# both Viking landing sites (years 1-2 for VL-1 and year 1 for VL-2).  This uses data
# from Neal Johnson, JGR 98 p. 10,973

amp=create(x=5,y=1,z=1,start=0,step=0,format=float)
phase=create(x=5,y=1,z=1,start=0,step=0,format=float)

amp[1,,]=0.706
amp[2,,]=0.581
amp[3,,]=0.108
amp[4,,]=0.061
amp[5,,]=0.016

phase[1,,]=92.27
phase[2,,]=-130.80
phase[3,,]=-69.76
phase[4,,]=-10.00
phase[5,,]=49.55

# Find relative pressure compared to mean

p_relative= amp * sin(2*3.1415927 * frac + phase * 3.1415927/180)

p_relative=sum(p_relative)+1


# Find absolute pressure using 6.1 mbar as mean surface pressure at 0 km.
# Use scale height of 10 km

scale_height=-($1)/10.

p_surf=6.1 * p_relative * exp(scale_height)

return(p_surf)

}








define read_vm(axis,style) {

if ($ARGC == 0) {
	printf (" \n")
	printf (" Read vm format text file. \n")
	printf (" $1 = filename \n")
	printf (" $2 = 0 - No xaxis (optional) \n")	
	printf (" axis = direction to store spectra (\"y\" or Default=\"x\" \n")
	printf (" style = Name of sample name in structure (1=\"label\" (Default) or 2=\"sample_name\"\n")
	printf (" \n")
	return(0)

}

filename=$1
direction="x"
label_style=1
if(HasValue(axis)) direction=axis
if(HasValue(style)) label_style=style

nvecs=atoi(syscall("less \"" + filename +"\" ")[7:,1])

out=struct(data,label,xaxis)

if ($ARGC == 1) {

	out.label=syscall("less \"" + filename+"\" ")[,3:(1+nvecs)]
	out.data=ascii(filename,format=float,row=(nvecs+1))
	out.xaxis=translate(out.data[1,,],from=y,to=z)
	out.data=translate(out.data[2:(nvecs),,],from=y,to=z)
}


if ($ARGC == 2) {

	out.label=syscall("less " + filename+" ")[,2:(1+nvecs)]
	out.data=ascii(filename,format=float,row=(nvecs+1))
	out.xaxis=0
	out.data=translate(out.data[1:(nvecs),,],from=y,to=z)
}

if(direction=="y"){
	out.data=translate(out.data,from=x,to=y)
}
if(label_style==2) {
	out.sample_name=out.label
	remove_struct(out,"label")
}


return(out)

}


define write_vm() {

if ($ARGC == 0) {
	printf (" \n")
	printf (" Write standard spectral structure to a vm text file. \n")
	printf (" $1 = structure \n")
	printf (" $2 = filename \n")
	printf (" \n")
	return(0)

}

out=$1

if (dim(out.data)[2,,]>1 && dim(out.xaxis)[3,,]!=1) {

        out.data=translate(out.data,from=x,to=y)
}

nvecs=dim(out.data)[1,,]
verbose=0
if (HasValue(out.sample_name)==1 && HasValue(out.label)==0) {

        out.label=out.sample_name
}
verbose=2

if (dim(out.xaxis)[3,,] > 1) {

	nvecs=nvecs+1
	out.label=cat("xscale",out.label,axis=y)
}

out.label=cat("label "+nvecs+"",out.label,axis=y)

write(out.label,"temp_label",type=ascii,force=1)

if (dim(out.xaxis)[3,,] > 1) {

	out.data=translate(out.xaxis,from=z,to=y)//translate(out.data,from=z,to=y)
}

write(out.data,"temp_data",type=ascii,force=1)
syscall("cat temp_label temp_data > "+$2+"")
syscall("rm temp_label temp_data")
return(0)

}



define make_mtesx() {
# make correct Mini-TES 
# usage: mtx = make_mtesx()  [-options]

# p. christensen after TES II make_xt2  6/22/02
# spectra compared to ASUlib quartz and dolomite on 4/9/03 to verify xaxis position
#  there had been some confusion and uncertainty about 1st point in spectrum returned
#  from FFT being DC term, but the algorithm below returns an xaxis that matches the ASU 
#  library.  PRC

npts = 512

start_pt = 35
end_pt = 201

# set wavenumber of first point downlinked (14 or 28 in PROM) #

# set prime factor
nprime = 1024;

laser_wavelength = 0.978
delta_x = (1./(nprime * laser_wavelength))*10000.
#printf(" delta_x = %.2f\n", delta_x)

mtx = create(npts, 1, 1, format=float)

for(i=1; i<= npts; i+=1) {
	# index (i) of 1 = DC term in FFT; index 2 = 9.9853 cm-1; etc
	# index 1 = channel 0; etc
	mtx[i] = float(i-1) * delta_x
}

# truncate after 4 decimal places 

mtx = mtx * 10000.
mtx = int(mtx)
mtx = mtx/10000.

mtx = translate(mtx,x,z)[,,start_pt:end_pt]
return(mtx)

}

define make_tesx(scan_len) {
# return tes 143-pt wavenumber values (or 286-pt) ("x-axis")
# p. christensen 4/29/97 
# double added prc 12/6/07

# copied from make_xt2 without change
# 4/19/07 p. christensen

if(HasValue(scan_len)==0) {
	scan_len = 1
}
if(scan_len == 2) {
	# double scan
	npts=286
	first_pt=28;
	# set prime factor to default (1344; Det. 2 & 5) */
	nprime = 1344;
	# tes resolution for 1344 pt fft
	tes_res = 10.5808820/2.
} else if(scan_len == 1) {
	# single scan
	npts=143
	first_pt=14;
	# set prime factor to default (1344; Det. 2 & 5) */
	# nprime = 1344;
	# tes resolution for 1344 pt fft
	tes_res = 10.5808820
} else {
	# bad value of scan len
	return("bad value of scan len")
}

xaxis=create(1,1,npts,start=tes_res*first_pt, step = tes_res, format=float)

/* truncate after 4 decimal places -9/4/97 */
xaxis=xaxis*10000
xaxis=int(xaxis)
xaxis=xaxis/10000.

return(xaxis)

}

define make_lab1x() {
#added $DV_EX support
#Added $DV_SCRIPT_FILES support - 5-13-08
#added instrument_parameters directory to $DV_SCRIPT_FILES 1-30-09

a = read($DV_SCRIPT_FILES+"/instrument_parameters/lab1_xaxis.vicar")
return(a)
}

define make_themisx() {
# Make 10 point THEMIS wavenumber scale
# modified from themisx; reverse order, put values in z not y
# themisx in /u/bandfield/josh.dvrc:3634
# p. christensen 4/19/07

xthemis=create(1,1,10, format=float)

xthemis[1,1,1]=675.3144
xthemis[1,1,2]=796.4805
xthemis[1,1,3]=850.8455
xthemis[1,1,4]=909.8835
xthemis[1,1,5]=983.6909
xthemis[1,1,6]=1075.0894
xthemis[1,1,7]=1174.8452
xthemis[1,1,8]=1267.6153
xthemis[1,1,9]=1477.8360
# band 10 value made slightly different from band 9 for plotting convenience
xthemis[1,1,10]=1477.85
#xthemis[1,1,10]=1477.8360
return(xthemis)
}

define lsubs() {

if ($ARGC == 0) {
        printf ("\n")
        printf ("Give approximation of Lsubs given et\n")
	return(0)
}

const=59340000.0
et=$1
pi=3.1415927
scale=et/const*2*pi+2

lsubs=(et+46090000)%const*360.0/const+sin(et*2*pi/const+0.4)*10.75+6*sin(-2*scale-0.3)/9.-pow(sin(scale)+1,4)/9-3*et/2000000000.0-8.75

return(lsubs)
}


define deconstruct() {

if ($ARGC == 0) {
	printf (" \n")
	printf (" Reduce an image into a single strip array: \n")
	printf (" $1 = input image \n")
	printf (" \n")
	return(0)

}


cube=$1

vecs=dim(cube)


em=translate(cube[,1,],from=x, to=y)
em=translate(em, from=z, to=x)

count=2

while (count <= vecs[2,,]) {

	em_temp=translate(cube[,(count),],from=x, to=y)
	em_temp=translate(em_temp,from=z, to=x)
	em=cat(em,em_temp,y)
	count=count+1
}

return(em)
}



define rm_quote() {
if ($ARGC == 0) {
	printf("Remove single or double quotes from an array of labels and replaces it with nothing\n")
	printf(" usage: a = remove_quote(label_array)\n")
	printf(" \n")
	printf(" where:\n")
	printf(" 	label_array is an array of labels (assumed to be in y-axis)\n")
	printf("\np. christensen 11/24/2010\n")
	return(0)
}

a = $1

b=strsub(a,"\'","")
b=strsub(b,"\"","")

return(b)

}



define bs() {

	#returns a backslash string nomatter what the input.
	return(string(byte(92)))
}