############################################################
### filename                src_rotation_example.lsf
### Update source values as required
### This script takes the fields as returned by the source in XYZ coordinates.
### It defines a coordinate transformation from this global system to the sources coordinates 
### as described by IA, phi, and theta. 
### This may be useful in working out your own coodinate transformation



### Run to ensure the fields are updated after any changes. May not be necessary
switchtolayout;
run;
selectpartial('source');

######################################
###Construct Polarization Matrix
######################################

pol = get('polarization angle');
?"Source polarization angle = " + num2str(pol);
pol = pol*(pi/180);

P = [cos(pol), -sin(pol), 0; sin(pol), cos(pol), 0; 0, 0, 1];

###Construct Elevation Rotation Matrix
theta = get('angle theta');
?"Source elevation angle = " + num2str(theta);
theta = theta*(pi/180);

A =  [cos(theta), 0, sin(theta); 0, 1, 0; -sin(theta), 0, cos(theta)];

###Construct Horizontal Rotation Matrix
phi = get('angle phi');
?"Source horizontal angle = " + num2str(phi);
phi = phi*(pi/180);

B = [cos(phi), -sin(phi), 0; sin(phi), cos(phi), 0; 0, 0, 1];

###Construct injection axis matrix
IA = get('injection axis');
if(IA == 'x-axis'){
    ?"Source injection axis is " + IA;
    C = [0, 0, 1; 1, 0, 0; 0, 1, 0];
} else if(IA == 'y-axis'){
    ?"Source injection axis is " + IA;
    C = [0, 1, 0; 0, 0, 1; 1, 0, 0];
} else if (IA == 'z-axis'){
    ?"Source injection axis is " + IA;
    C = [1, 0, 0; 0, 1, 0; 0, 0, 1];
}

### Consider direction of propogation
dirc = get('direction');
if(dirc == 'Backward'){C = -1*C;} #This is a reflection about the origin

### Final Matrix Rotation from global to local coordinates
#D = mult(transpose(P),mult(transpose(A),transpose(B),transpose(C))); # This transformation will leave E X polarized in E_local
D = mult(transpose(A),transpose(B),transpose(C)); ##Exclude P if want the fields rotated in the plane by pol angle

#####################################
### Get source fields
#####################################

SF = getresult('source', 'fields');
SF_E=pinch(SF.E);
x_G = pinch(SF.x);
y_G = pinch(SF.y);
z_G = pinch(SF.z);

X = meshgridx(x_G,y_G);
Y = meshgridy(x_G,y_G);
##Init matrices
xyz = XYZ = E_local = E_global = matrix(length(SF_E(:,1,1)),length(SF_E(1,:,1)),length(SF_E(1,1,:)));

###These are the E fields in global coordinate system
E_global(:,:,1) = SF_E(:,:,1);
E_global(:,:,2) = SF_E(:,:,2);
E_global(:,:,3) = SF_E(:,:,3);

#mesh grid of source plane in local coordinates
XYZ(:,:,:) = [X,Y,matrix(length(x_G),length(y_G))];

### Transform global to local coordinates using D
for (k=1:length(pinch(E_local(:,1,1)))){
E_local(k,:,:) = transpose(mult(D,transpose(pinch(E_global(k,:,:)))));
xyz(k,:,:) = transpose(mult(D,transpose(pinch(XYZ(k,:,:))))); #This is the plane of constant phase
}

########################################
### Reconstruct datasets
########################################
E_XYZ = rectilineardataset("E",x_G,y_G,z_G);
E_XYZ.addattribute("E global",E_global(:,:,1),E_global(:,:,2),E_global(:,:,3));

E_xyz = rectilineardataset("E",x_G,y_G,z_G);
E_xyz.addattribute("E local",E_local(:,:,1),E_local(:,:,2),E_local(:,:,3));

### Visualize
visualize(D); ##Rotation matrix
visualize(E_XYZ); ##Fields in the global coordinates
visualize(E_xyz); ##Field in coordinate system with k||z