% geo2mag3.m is used to convert geographical coords to
% magnetic coords. It uses the values given in Appendix 3,
% p534, of Introduction to Space Physics. NOTE that input
% and output are in rads, not degrees. Unlike geo2mag.m and
% geo2mag2.m, it allows any size input matrix and returns 
% the same size output.
%%
%% Author: Phillip A. Webb
%% Last updated: June 9 2004                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           
%%

function [lonmg,latmg]=geo2mag3(longo,latgo)

[yi,xi]=size(longo);                                      %size of input matrices      
longo=longo(:);                                           %convert longo to column vector
latgo=latgo(:);                                           %convert latgo to column vector

radgo=ones(size(longo));                                  %unit radi vector
[xv,yv,zv]=sph2cart(longo,latgo,radgo);                   %convert from spherical to cartesian coords
tm1=[0.32110 -0.92756 -0.19112];
tm2=[0.94498 0.32713 0];
tm3=[0.06252 -0.18060 0.98157];
tm=[tm1;tm2;tm3];                                         %transformation matrix
res=tm*[xv';yv';zv'];                                     %transformed result matrix
[lonmg,latmg,radmg]=cart2sph(res(1,:),res(2,:),res(3,:)); %convert from cartesian to spherical

lonmg=reshape(lonmg,yi,xi);                               %required longitude values
latmg=reshape(latmg,yi,xi);                               %required latitude values

