function cktiles(shot,t)
 
 if nargin == 0
  cktiles(47650,1.0);
  cktiles(40624,0.8);
  cktiles(47651,1.0);
  cktiles(45109,1.0);
  return
 end
 
 % tile contour
%   [xt,yt] = tilec([1 2],0); nt = 32;
[xt,yt] = tilec([1 2 3 4],2); nt = 16;
xt = xt(:)'; yt = yt(:)';
 rt = sqrt(xt.^2+yt.^2); at = atan2(yt,xt);
 
 r0 = min(rt(:)); z0 = 0.5494 * linspace(-1,1,1001); ls = 0.003; Ps = 1e6;
 
 mdsopen(shot)
 rbt = mdsdata('\MAGNETICS::RBPHI[$]',t);
 psi = psitbxtcv(shot,t,'+0');
 fsd = psitbxp2p(psi,'FS');
 
 % power flux and flux expansion at midplane separatrix
 rs = fsd.rmag + fsd.x(end,1);
 qs = Ps/(4*pi*rs*ls);
 drpsis = psitbxf2f(psi,psitbxgrid('C','P',{rs fsd.zmag}),[1 0]); drpsis = drpsis.x;
 
 % find vertical position of strike point on reference cylinder
 gr = psitbxgrid('C','G',{r0 z0});
 psir = psitbxf2f(psi,gr); psir =   psir.x;
 k = find(psir(1:end-1).*psir(2:end)<0,1);
 zd = z0(k) - (z0(k+1)-z0(k))/(psir(k+1)-psir(k))*psir(k);
 
 % power flux and magnetic field at strike point on reference cylinder
 gd = psitbxgrid('C','G',{r0 zd});
 drpsid = psitbxf2f(psi,gd,[1 0]); drpsid = drpsid.x;
 dzpsid = psitbxf2f(psi,gd,[0 1]); dzpsid = dzpsid.x;
 qd = qs*rs/r0*dzpsid/drpsis;
 ld = ls*drpsis/dzpsid;
 brd = -dzpsid/(2*pi*r0);
 bzd =  drpsid/(2*pi*r0);
 btd = abs(rbt)/r0;
 
 % projection and shadow of tile on reference cylinder
 ar = at -    btd/brd*log(rt/r0);
 dar = diff(ar);
 kr = true(size(ar));
 for k1 = find(dar(1:end-1)>0 & dar(2:end)<0)+1
   kr(k1:end) = ar(k1:end) > ar(k1);
 end
 % only one period
 xt = xt(end/2+1:end);
 yt = yt(end/2+1:end);
 rt = rt(end/2+1:end);
 at = at(end/2+1:end);
 ar = ar(end/2+1:end);
 kr = kr(end/2+1:end);
 dt = sqrt(diff(xt).^2+diff(yt).^2);
 lt = [0 cumsum(dt)];
 zt = linspace(-1,7,length(rt))'*ld; dzt = mean(diff(zt));
 zr = zt-r0*bzd/brd*log(rt/r0);
 
 % power flux projection
 qt = -qd*ld*(diff(exp(-max(zr,0)/ld)));
 qt = (qt(:,1:end-1)+qt(:,2:end))/2.*(r0*diff(ar));
 qt = qt./(dt*dzt);
 Pt = nt*sum(sum(qt.*(dt*dzt)));
 qt(end+1,:) = 0; qt(:,end+1) = 0; % for pcolor
 qt(:,~kr) = 0;
 
 % plot
 tag = sprintf('#%d %4.0s',shot,t);
 h = findobj('tag',tag);
 if h, figure(h)
 else figure('tag',tag)
 end
 clf
 subplot(222)
 set(pcolor(lt,zt,qt/1e6),'edgecolor','none'), axis equal, xlim(lt([1 end])), colorbar
 hold on, plot(lt,max(qt)/max(qt(:))*zt(end),'-k',lt,rt-rt(1),'-k'), hold off
 title(sprintf('z_d=%5.2fm, \\lambda_d=%3.1fmm - P_t=%3.0fkW - q_t<%2.0fMW/m^2',zd,ld*1e3,Pt/1e3,max(qt(:))/1e6))
 subplot(221)
 title(sprintf('#%d %4.2fs - P_s=%1.0fMW, \\lambda_s=%1.0fmm',shot,t,Ps/1e6,ls*1e3))
 hold on, tcvview('tvp',shot,t), hold off
 colormap(flipud(hot))
end

function [xt,yt] = tilec(ind,port)
 r1 = 15e-3; r2 = 40e-3; r3 = 2e-3; l1 = 72.2e-3-r3; l2 = 4.4e-3;
 xr = 3e-3; yr = 50.8e-3; ar = 3.07/180*pi; rr = 1154e-3; rv = 1160e-3;
 a2 = -2.561/180*pi; a1 = 360/32/180*pi;
 n = [100 130 100 7 200];
 
 % tile coordinates
 t1 = linspace(-pi/2,0,n(1));
 t2 = linspace(0,pi/2,n(3));
 xt = [   r1*cos(t1) repmat(r1,1,n(2)-2)       r1-r3+r3*cos(t2) r1-r3-(1:n(4))/n(4)*l2];
 yt = [r2+r2*sin(t1) (1:n(2)-2)*l1/(n(2)-1)+r2 r2+l1+r3*sin(t2) repmat(r2+l1+r3,1,n(4))];
 
  % rail coordinates
 tmp = [cos(ar) sin(ar);-sin(ar) cos(ar)]*[xt;yt];
 xt = tmp(1,:); yt = tmp(2,:);
 xt = xt-xr; yt = yt-yr;
 
 % TCV cartesian coordinates
 [xt,yt] = rot(rr-xt,yt,a2,-1);
 
 % add port
 switch port
  case 0
  case 1
   a = (n(5):-1:0)*(a1/(2*n(5)+1));
   xt = [rv*cos(a) xt];
   yt = [rv*sin(a) yt];
  case 2
   aa = [a1*n(5)/(2*n(5)+1) a1/4 a1/4 0];
   rr = [rv+24e-3 rv+24e-3 rv rv];
   a = []; r = [];
   for k = 1:length(aa)-1
    tmp = linspace(aa(k),aa(k+1),n(5)/2);
    a = [a tmp(1:end-1)];
    tmp = linspace(rr(k),rr(k+1),n(5)/2);
    r = [r tmp(1:end-1)];
   end
   a = [a aa(end)];
   r = [r rr(end)];
   xt = [r.*cos(a) xt];
   yt = [r.*sin(a) yt];
   
 end
 
 % rotate and duplicate
 ind = ind(:)';
 yf = 1-2*mod(ind,2);
 a = a1*floor(ind/2);
 if port, a = a1*ind; end
 [xt,yt] = rot(xt(:),yt(:),a,yf);
 k = ~logical(mod(ind,2));
 xt(:,k) = flipud(xt(:,k));
 yt(:,k) = flipud(yt(:,k));
 
 function [X,Y] = rot(x,y,a,yf)
  X = x*cos(a) - y*(yf.*sin(a));
  Y = x*sin(a) + y*(yf.*cos(a));
 end
end