function md=radarpower(md,varargin)
%RADARPOWER - overlay a power radar image on an existing mesh
%
%   This routine will overlay a power radar image on an existing mesh.
%   The power amplitude will be output to vel for now.
%   In the future, think about a field to hold this value.
%
%   Usage:
%      md=radarpower(md,options);
%      md=radarpower(md)

%If gdal does not work, uncomment the following line
%setenv('LD_LIBRARY_PATH','/proj/ice/larour/issm/trunk/externalpackages/gdal/install/lib/');
%Parse inputs
if nargin==1,
	options=pairoptions;
else
	options=varargin{:};
	if ~isa(options,'pairoptions'),
		options=pairoptions(varargin{:});
	end
end

highres=getfieldvalue(options,'highres',0);
xlim=getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]);
ylim=getfieldvalue(options,'ylim',[min(md.mesh.y) max(md.mesh.y)]);
posting=getfieldvalue(options,'posting',0); % 0 -> image posting default

%find gdal coordinates
x0=min(xlim); x1=max(xlim);
y0=min(ylim); y1=max(ylim);

%figure out if we should go look for Greenland or Antarctica geotiff, or if user provided one.
if ~exist(options,'overlay_image'),
	if strcmpi(md.mesh.hemisphere,'n'),
		if ~exist([issmdir() '/projects/ModelData/MOG/mog150_greenland_map.jpg']),
			error(['radarpower error message: file ' issmdir() '/projects/ModelData/MOG/mog150_greenland_map.jpg not found.']);
		end
		jpgim=[issmdir() '/projects/ModelData/MOG/mog150_greenland_map.jpg'];
		geom=load([issmdir() '/projects/ModelData/MOG/mog150_greenland_map.jpgw'],'ascii');

		%geom:   xposting nbcols nbrows yposting xmin ymax
		xmin=max(geom(5),x0);
		xmax=min(geom(5)+geom(1)*geom(2),x1);
		ymin=max(geom(6)-geom(4)*geom(3),y0);
		ymax=min(geom(6),y1);

		firstcol=max(1,floor((xmin-geom(5))/geom(1))); %x min
		firstrow=max(1,floor((geom(6)-ymax)/geom(4))); %y max
		numcols=floor((xmax-xmin)/geom(1)); % x posting
		numrows=floor((ymax-ymin)/geom(4)); % y posting
		pixelskip=max(1,ceil(posting/geom(1)));

		%Read and crop file
		disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
		im=imread(jpgim);
		im=im(firstrow:firstrow+numrows-1,firstcol:firstcol+numcols-1);
		md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
		md.radaroverlay.x=(xmin:(xmax-xmin)/(size(md.radaroverlay.pwr,2)-1):xmax);
		md.radaroverlay.y=(ymin:(ymax-ymin)/(size(md.radaroverlay.pwr,1)-1):ymax);

	elseif strcmpi(md.mesh.hemisphere,'s'),
		if highres,
			if ~exist([issmdir() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']),
				error(['radarpower error message: file ' issmdir() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif not found.']);
			end
			geotiff_name=[issmdir() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif'];
		else
			if ~exist([issmdir() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']),
				error(['radarpower error message: file ' issmdir() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif not found.']);
			end
			geotiff_name=[issmdir() '/projects/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif'];
		end

		%Name of image
		inputname='./temp.tif';
		eval(['!gdal_translate -quiet -projwin ' num2str(x0) ' ' num2str(y1) ' ' num2str(x1) ' ' num2str(y0) ' ' geotiff_name ' ' inputname ]);

		%Read in temp.tif:
		im=imread('temp.tif','TIFF');
		pixelskip=max(1,ceil(posting/((x1-x0)/(size(im,2)))));
		md.radaroverlay.pwr=double(flipud(im(1:pixelskip:end,1:pixelskip:end)));
		md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
		md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);

		%Erase image
		system('rm -rf ./temp.tif');

	else
		error('field hemisphere should either be ''n'' or ''s''');
	end
else
	%ok, user provided an image. check we also have overlay_xlim and overlay_ylim  options, to know what range of coordinates the image covers.
	if (~exist(options,'overlay_xlim') | ~exist(options,'overlay_xlim')| ~exist(options,'overlay_xposting')| ~exist(options,'overlay_yposting')),
		error('radarpower error message: please provide overlay_xlim, overlay_ylim, overlay_xposting and overlay_yposting options together with overlay_image option');
	end
	overlay_image=getfieldvalue(options,'overlay_image');
	overlay_xlim=getfieldvalue(options,'overlay_xlim');
	overlay_ylim=getfieldvalue(options,'overlay_ylim');
	overlay_xposting=getfieldvalue(options,'overlay_xposting');
	overlay_yposting=getfieldvalue(options,'overlay_yposting');

	sizex=floor((x1-x0)/overlay_xposting);
	sizey=floor((y1-y0)/overlay_yposting);
	topleftx=floor((x0-overlay_xlim(1))/overlay_xposting); % x min
	toplefty=floor((overlay_ylim(2)-y1)/overlay_yposting); % y max


	%Read and crop file
	disp('Warning: expecting coordinates in polar stereographic (Std Latitude: 70ºN Meridian: 45º)');
	im=imread(overlay_image);
	im=im(toplefty:toplefty+sizey,topleftx:topleftx+sizex);
	md.radaroverlay.pwr=double(flipud(im));
	md.radaroverlay.x=(x0:(x1-x0)/(size(md.radaroverlay.pwr,2)-1):x1);
	md.radaroverlay.y=(y0:(y1-y0)/(size(md.radaroverlay.pwr,1)-1):y1);
end
