%
% This is a copy of the program pgdemo4 from the pgplot distribution,
% translated into octave/matlab notation.  I translated it pretty
% literally; I don't really know what most of it does.
% Some of this could be more compactly and elegantly expressed in octave if I
% bothered to try.
%
1;

%-----------------------------------------------------------------------
% Set the viewport, allowing margins around the edge for annotation.
% (This is similar in effect to PGVSTD, but has different margins.)
% The routine determines the view-surface size and allocates margins
% as fractions of the minimum of width and height.
%-----------------------------------------------------------------------
function setvp
  pgsvp(0, 1, 0, 1);
  [vpx1, vpx2, vpy1, vpy2] = pgqvp(1);
  d = min(vpx2-vpx1, vpy2-vpy1)/40;
  vpx1 = vpx1 + 5*d;
  vpx2 = vpx2 - 2*d;
  vpy1 = vpy1 + 8*d;
  vpy2 = vpy2 - 2*d;
  pgvsiz(vpx1, vpx2, vpy1, vpy2);
endfunction

%-----------------------------------------------------------------------
% Draw the enclosing rectangle of the subarray to be contoured,
% applying the transformation TR.
%
% For a contour map, the corners are (I1,J1) and (I2,J2); for
% a gray-scale map, they are (I1-0.5,J1-0.5), (I2+0.5, J2+0.5).
%-----------------------------------------------------------------------
function outlin(i1,i2,j1,j2,tr)
  xw = [i1; i1; i2; i2; i1];
  yw = [j1; j2; j2; j1; j1];
  for k=1:5
    t = xw(k);
    xw(k) = tr(1) + tr(2)*t + tr(3)*yw(k);
    yw(k) = tr(4) + tr(5)*t + tr(6)*yw(k);
  end

  pgline(xw, yw);
endfunction

%-----------------------------------------------------------------------
% Set a "palette" of colors in the range of color indices used by
% PGIMAG.
%-----------------------------------------------------------------------
function palett(type)
  GL = [0; 1];
  GR = [0; 1];
  GG = [0; 1];
  GB = [0; 1];

  RL = [-0.5; 0.0; 0.17; 0.33; 0.50; 0.67; 0.83; 1.0; 1.7];
  RR = [ 0.0; 0.0;  0.0;  0.0;  0.6;  1.0;  1.0; 1.0; 1.0];
  RG = [ 0.0; 0.0;  0.0;  1.0;  1.0;  1.0;  0.6; 0.0; 1.0];
  RB = [ 0.0; 0.3;  0.8;  1.0;  0.3;  0.0;  0.0; 0.0; 1.0];

  HL = [0.0; 0.2; 0.4; 0.6; 1.0];
  HR = [0.0; 0.5; 1.0; 1.0; 1.0];
  HG = [0.0; 0.0; 0.5; 1.0; 1.0];
  HB = [0.0; 0.0; 0.0; 0.3; 1.0];

  WL = [0.0; 0.5; 0.5; 0.7; 0.7; 0.85; 0.85; 0.95; 0.95; 1.0];
  WR = [0.0; 1.0; 0.0; 0.0; 0.3;  0.8;  0.3;  1.0;  1.0; 1.0];
  WG = [0.0; 0.5; 0.4; 1.0; 0.0;  0.0;  0.2;  0.7;  1.0; 1.0];
  WB = [0.0; 0.0; 0.0; 0.0; 0.4;  1.0;  0.0;  0.0; 0.95; 1.0];

  AL = [0.0; 0.1; 0.1; 0.2; 0.2; 0.3; 0.3; 0.4; 0.4; 0.5; 0.5; ...
	0.6; 0.6; 0.7; 0.7; 0.8; 0.8; 0.9; 0.9; 1.0];
  AR = [0.0; 0.0; 0.3; 0.3; 0.5; 0.5; 0.0; 0.0; 0.0; 0.0; 0.0; ...
	0.0; 0.0; 0.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0];
  AG = [0.0; 0.0; 0.3; 0.3; 0.0; 0.0; 0.0; 0.0; 0.8; 0.8; 0.6; ...
	0.6; 1.0; 1.0; 1.0; 1.0; 0.8; 0.8; 0.0; 0.0];
  AB = [0.0; 0.0; 0.3; 0.3; 0.7; 0.7; 0.7; 0.7; 0.9; 0.9; 0.0; ...
	0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0];

  TL = [0.0; 0.5; 0.5; 1.0];
  TR = [0.2; 0.6; 0.6; 1.0];
  TG = [0.0; 0.0; 0.5; 1.0];
  TB = [1.0; 0.0; 0.0; 0.0];

  switch (type)
    case (1)			%        -- gray scale
         pgctab(GL, GR, GG, GB, 1.0, 0.5);
    case (2)			%        -- rainbow
         pgctab(RL, RR, RG, RB, 1.0, 0.5);
    case (3)			%        -- heat
         pgctab(HL, HR, HG, HB, 1.0, 0.5);
    case (4)			%        -- weird IRAF
	 pgctab(WL, WR, WG, WB, 1.0, 0.5);
    case (5)			%        -- AIPS
         pgctag(AL, AR, AG, AB, 1.0, 0.5);
    case (6)			%        -- TJP
         pgctab(TL, TR, TG, TB, 1.0, 0.5)
  endswitch
endfunction

mxi = 40;
mxj = 40;
i = 1:mxi; j = 1:mxj;
[jm, im] = meshgrid(i, j);		% Note order reversal!
f = cos(0.6*sqrt(im*80/mxi)-16*jm/(3*mxj)).*cos(16*im/(3*mxi))+(im/mxi-jm/mxj);

fmax = max(max(f));
fmin = min(min(f));
angle = -120*pi/180;		% Image rotation angle in radians.
radius = 40*sqrt(2);
tr(1) = -(40/mxi)*cos(angle)+(40/mxj)*sin(angle) - radius*cos(angle+45*pi/180);
tr(2) = 80*cos(angle)/mxi;
tr(3) = -80*sin(angle)/mxj;
tr(4) = -(40/mxi)*sin(angle)-(40/mxj)*cos(angle) - radius*sin(angle+45*pi/180);
tr(5) = 80*sin(angle)/mxi;
tr(6) = 80*cos(angle)/mxj;

%
% Clear the screen; set up window and viewport.
%
pgpage;
setvp;
pgwnad(-40, 40, -40, 40);
pgsci(1);

%
% Draw the map with PGGRAY.  
%
pggray(f,1,mxi,1,mxj,fmax,fmin,tr)

%
% Overlay contours in red:
%
pgsci(2);

for i = 1:21
  alev = fmin + (i-1)*(fmax-fmin)/20;
  if (rem(i, 5) == 0)
    pgslw(3);
  else
    pgslw(1);
  end

  if (i < 10)
    pgsls(2);
  else
    pgsls(1);
  end
  pgcont(f, 1, mxi, 1, mxj, alev, -1, tr);
end
pgsls(1);
pgslw(1);

%
% Annotate the plot:
%
pgsci(2);
outlin(1, mxi, 1, mxj, tr);
pgsci(5);
pgmtxt('t', 1, 0, 0, 'Routimes PGGRAY, PGCONT, PGWEDG');
pgbox('bcnts', 0, 0, 'bcnts', 0, 0);

%
% Draw a wedge.
%
pgsch(0.8);
pgwedg('BG', 3, 4, fmax, fmin, 'Elevation');
pgsch(1);

%-----------------------------------------------------------------------
% PGIMAG
%-----------------------------------------------------------------------
% Clear the screen. Set up window and viewport.
%
pgpage;
setvp;
pgwnad(-40, 40, -40, 40);
pgsci(1);
%
% Set up the color map.
%
palett(2);

%
% Draw the map with PGIMAG.  
%
pgimag(f,1,mxi,1,mxj,fmin,fmax,tr);
%
% Overlay contours in white.
%
pgsci(1)
for i=1:21
  alev = fmin + (i-1)*(fmax-fmin)/20.0;
  if (rem(i,5) == 0)
      pgslw(3);
  else
      pgslw(1);
  end
  if (i < 10)
      pgsls(2);
  else
      pgsls(1);
  end
  pgcont(f,1,mxi,1,mxj,alev,-1,tr);
end

pgsls(1);
pgslw(1);
%
% Annotate the plot.
%
pgsci(2);
outlin(1,mxi,1,mxj,tr);
pgsci(5);
pgmtxt('t',1.0,0.0,0.0,'Routines PGIMAG, PGCONT, PGWEDG');
pgbox('bcnts',0.0,0,'bcnts',0.0,0);
%
% Draw a wedge.
%
pgsch(0.8);
pgwedg('BI', 3.0, 4.0, fmin, fmax, 'Elevation');
pgsch(1.0);