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.6 2005/04/11 18:15:18 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 % load External
0024 % extlat = (pi/2)-acos(external(:,1));
0025 % extlong = external(:,2);
0026 
0027 % Observatory data
0028 
0029 CordntH = getifo('H1');
0030 
0031 lhphi = CordntH.LAT;
0032 lhpsi = CordntH.LONG; 
0033 
0034 lhx = cos(lhpsi)*cos(lhphi); 
0035 lhy = sin(lhpsi)*cos(lhphi);
0036 lhz = sin(lhphi);
0037 
0038 CordntL = getifo('L1');
0039 
0040 llphi = CordntL.LAT; 
0041 llpsi = CordntL.LONG; 
0042 
0043 llx = cos(llpsi)*cos(llphi); 
0044 lly = sin(llpsi)*cos(llphi);
0045 llz = sin(llphi);
0046 
0047 CordntL = getifo('TAMA');
0048 
0049 tmphi = CordntL.LAT; 
0050 tmpsi = CordntL.LONG; 
0051 
0052 tmx = cos(llpsi)*cos(llphi); 
0053 tmy = sin(llpsi)*cos(llphi);
0054 tmz = sin(llphi);
0055 
0056 CordntL = getifo('GEO');
0057 
0058 gophi = CordntL.LAT; 
0059 gopsi = CordntL.LONG; 
0060 
0061 gox = cos(llpsi)*cos(llphi); 
0062 goy = sin(llpsi)*cos(llphi);
0063 goz = sin(llphi);
0064 
0065 CordntL = getifo('VIRGO');
0066 
0067 vgphi = CordntL.LAT; 
0068 vgpsi = CordntL.LONG; 
0069 
0070 vgx = cos(llpsi)*cos(llphi); 
0071 vgy = sin(llpsi)*cos(llphi);
0072 vgz = sin(llphi);
0073 
0074 % Map data
0075 zoomfac = 1.4;
0076 
0077 load coast;
0078 
0079 x = [];
0080 y = [];
0081 z = [];
0082 
0083 for i = 1:length(long)
0084    phi(i) = lat(i)*pi/180;
0085    psi(i) = long(i)*pi/180;
0086    
0087    x(i) = cos(psi(i))*cos(phi(i)); 
0088    y(i) = sin(psi(i))*cos(phi(i));
0089    z(i) = sin(phi(i));
0090 end
0091 
0092 Llo.globe = [llx; lly; llz];
0093 Llo.map = [llpsi; llphi];
0094 Lho.globe = [lhx; lhy; lhz];
0095 Lho.map = [lhpsi; lhphi];
0096 Coast.globe = [x; y; z];
0097 Coast.map = [psi; phi];
0098 
0099 % Plotting
0100 % Globe
0101 figure(1)
0102 
0103 clf
0104 sphere(36);
0105 h1 = findobj('Type', 'surface');
0106 set(h1, 'CData', flipud(zz), 'FaceColor', 'texturemap')
0107 set(gca, 'CameraPosition', zoomfac*[lhx, lhy, lhz], 'CameraTarget', ...
0108     [0, 0, 0])
0109 axis off;
0110 axis square;
0111 colorbar;
0112 hold on;
0113 
0114 plot3(x(1:end), y(1:end), z(1:end), 'w');
0115 plot3(llx, lly, llz, 'kx');
0116 plot3(lhx, lhy, lhz, 'kx');
0117 plot3(gox, goy, goz, 'kx');
0118 plot3(tmx, tmy, tmz, 'kx');
0119 plot3(vgx, vgy, vgz, 'kx');
0120 
0121 % Map
0122 figure(2)
0123 
0124 clf
0125 pcolor(xx, yy, zz);
0126 shading('interp');
0127 colorbar('horiz');
0128 hold on;
0129 
0130 plot(psi(1:end), phi(1:end), 'w');
0131 %plot(extlong(1:end), extlat(1:end), '.k')
0132 plot(llpsi, llphi, 'kx');
0133 plot(lhpsi, lhphi, 'kx');
0134 plot(gopsi, gophi, 'kx');
0135 plot(tmpsi, tmphi, 'kx');
0136 plot(vgpsi, vgphi, 'kx');
0137 
0138 axis equal;
0139 axis([-pi, pi, -pi/2, pi/2]);
0140 xlabel('Longitude');
0141 ylabel('Latitude');
0142 xlim([-pi, pi]);
0143 ylim([-pi/2, pi/2]);
0144 set(gca, 'XTick', -pi:pi/4:pi);
0145 set(gca, 'XTickLabel', {'180 W', '135 W', '90 W', '45 W', '0', '45 E', ...
0146         '90 E', '135 E', '180 E'});
0147 set(gca, 'YTick', -pi/2:pi/4:pi/2);
0148 set(gca, 'YTickLabel', {'90 S', '45 S', '0', '45 N', '90 N'});
0149 shg
0150 
0151 return