ligomap

 LIGOMAP - utility to plot an array on world map with LIGO observatories
           marked

 ligomap(xx, yy, zz)

 xx    vector defining row dimension of array
 yy    vector defining column dimension of array
 zz    array to be ploted on world map (must be xx by yy in size)

 N.B. Needs data set coast.mat for world map

0001 function [Coast, Llo, Lho] = ligomap(xx, yy, zz)
0002 
0003 % LIGOMAP - utility to plot an array on world map with LIGO observatories
0004 %           marked
0005 %
0006 % ligomap(xx, yy, zz)
0007 %
0008 % xx    vector defining row dimension of array
0009 % yy    vector defining column dimension of array
0010 % zz    array to be ploted on world map (must be xx by yy in size)
0011 %
0012 % N.B. Needs data set coast.mat for world map
0013 
0014 % $Id: ligomap.m,v 1.4 2004/04/28 21:23:28 stuver Exp $
0015 
0016 % Special thanks to Derek Bridges for originally developing the
0017 % underpinnings of this function!
0018 
0019 Coast = deal(struct('globe', [], 'map', []));
0020 Llo = deal(struct('globe', [], 'map', []));
0021 Lho = deal(struct('globe', [], 'map', []));
0022 
0023 % Observatory data
0024 
0025 CordntH = getifo('H1');
0026 
0027 lhphi = CordntH.LAT;
0028 lhpsi = CordntH.LONG; 
0029 
0030 lhx = cos(lhpsi)*cos(lhphi); 
0031 lhy = sin(lhpsi)*cos(lhphi);
0032 lhz = sin(lhphi);
0033 
0034 CordntL = getifo('L1');
0035 
0036 llphi = CordntL.LAT; 
0037 llpsi = CordntL.LONG; 
0038 
0039 llx = cos(llpsi)*cos(llphi); 
0040 lly = sin(llpsi)*cos(llphi);
0041 llz = sin(llphi);
0042 
0043 % Map data
0044 zoomfac = 1.4;
0045 
0046 load coast;
0047 
0048 x = [];
0049 y = [];
0050 z = [];
0051 
0052 for i = 1:length(long)
0053    phi(i) = lat(i)*pi/180;
0054    psi(i) = long(i)*pi/180;
0055    
0056    x(i) = cos(psi(i))*cos(phi(i)); 
0057    y(i) = sin(psi(i))*cos(phi(i));
0058    z(i) = sin(phi(i));
0059 end
0060 
0061 Llo.globe = [llx; lly; llz];
0062 Llo.map = [llpsi; llphi];
0063 Lho.globe = [lhx; lhy; lhz];
0064 Lho.map = [lhpsi; lhphi];
0065 Coast.globe = [x; y; z];
0066 Coast.map = [psi; phi];
0067 
0068 % Plotting
0069 % Globe
0070 figure(1)
0071 
0072 clf
0073 sphere(36);
0074 h1 = findobj('Type', 'surface');
0075 set(h1, 'CData', flipud(zz), 'FaceColor', 'texturemap')
0076 set(gca, 'CameraPosition', zoomfac*[lhx, lhy, lhz], 'CameraTarget', ...
0077     [0, 0, 0])
0078 axis off;
0079 axis square;
0080 colorbar;
0081 hold on;
0082 
0083 plot3(x(1:end), y(1:end), z(1:end), 'w');
0084 plot3(llx, lly, llz, 'gx');
0085 plot3(lhx, lhy, lhz, 'wx');
0086 
0087 % Map
0088 figure(2)
0089 
0090 clf
0091 pcolor(xx, yy, zz);
0092 shading('interp');
0093 colorbar('horiz');
0094 hold on;
0095 
0096 plot(psi(1:end), phi(1:end), 'w');
0097 plot(llpsi, llphi, 'gx');
0098 plot(lhpsi, lhphi, 'wx');
0099 
0100 axis equal;
0101 axis([-pi, pi, -pi/2, pi/2]);
0102 xlabel('Longitude');
0103 ylabel('Latitude');
0104 xlim([-pi, pi]);
0105 ylim([-pi/2, pi/2]);
0106 set(gca, 'XTick', -pi:pi/4:pi);
0107 set(gca, 'XTickLabel', {'180 W', '135 W', '90 W', '45 W', '0', '45 E', ...
0108         '90 E', '135 E', '180 E'});
0109 set(gca, 'YTick', -pi/2:pi/4:pi/2);
0110 set(gca, 'YTickLabel', {'90 S', '45 S', '0', '45 N', '90 N'});
0111 shg
0112 
0113 return