% filename: current_center_deviations_v2b.m
% Script for Calculating Modular Coil Current Center Deviations.

%{
Created by Mark Smith
Princeton Plasma Physics Laboratory 
James Forrestal Campus 
P.O. Box 451 
Princeton, NJ 08543-0451 

Revised 1/07/08
..........................................................................
Validate Winding Coil Deviations Notes:
Two to three script files form the overall algorithm to create one excel file 
from the multiple Roamer Arm PowerInspect reports developed during the winding 
coil process.After this file is created, it is used to re-calculate the current 
center deviations.
This script file is used after locate_pwinspec_v.m, compile_pwinsp_v.m.
and possibly after additional manual data merging.
...........................................................................
Algorithm current_center_deviations_v2b.m notes:
The data/results reflect a ProE default coordinate system. 
The current center deviations dx (lateral to the septum) and dy (radial to the septum)
are determined at each clamp location. The deviations are calculated via averaging.
For example, each base data set consists of 4 measurments which are averaged
yielding 1 base value for this clamp location. After averaging at each clamp
location, there are 8 values representing the coil geometry (base,septum,top,side)
each for measurments sdies A and B. The current center is the center of 
Area for this cross-section which is computed via averaging of the 8 geometry
values. NOTE: side A septum and side values must be multipled -1 to account 
for sign convention. 
The current center deviations computed are then compared to Art Brooks' results. 
..........................................................................
Script notes: current_center_deviations_v2b.m
This code is used for B coils only.

INPUTS:
winding_data_B.xls 
The Excel tab labels and sequence are: baseA,baseB, septumA, septumB,topA,
topB, sideA, sideB.
The columns correspond to the following information.
column A:    clamp labels AND/OR sequence
columns B-D: X,Y,Z coordinate data
columns E-G: dx,dy,dz
column H:    3D deviation

OUTPUTS:
Comparison Plots.
Units are in inches.
%}

clc
clear

% Capture Current Matlab Directory.
working_matlab_dir=pwd;

[s1 s2]=uigetfile('*.xls','Select winding_data_B?.xls file.');

% Start program execution timer.
tic;

% Open matlab .mat file to obtain string variables for Excel tab labels.
tabs = struct2cell(open('Icenter_tabs.mat'));

% initialize variables
bsA=zeros(96,1);
bsB=bsA;
spA=bsA;
spB=bsA;
tpA=bsA;
tpB=bsA;
sdA=bsA;
sdB=bsA;
clamp_er=zeros(96,8);     % Track misssing/bad group/block data

% ************************************************************************
% Section (1): Process Excel Data file.
% Calculate averages at each clamp location for each Excel Sheet/Tab.
% ************************************************************************ 

for tbcc=1:8     % Excel tabs loop. Controls switch case.
    sheet=tabs{1}(1,tbcc);      % Excel Data labels.
    colA = zeros(3300,1);       % Initializing.
    colH = colA;                % Initializing.
    colA = xlsread(strcat(s2,s1),sheet{1},'A1:A3300'); % Read data col A.
    colH = xlsread(strcat(s2,s1),sheet{1},'H1:H3300'); % Read data col H.
    % Loop to search for each Clamp/Block #
    for blcc=1:96
       % Loop: search excel cells column A for block #.
       flg1=0;   % flag: check for bad excel group name
       for dtcc=1:length(colA)  % establishes row indice of excel.
           flg1=flg1+1;        
           if colA(dtcc)==blcc % Check if correct group/block found.
               flg1=0;                 % reset flag since group found.
               switch sheet{1}
                   % Case Bases/Tops read 4 excel rows of data
                   case 'baseA'
                       bsA(blcc) = mean(colH(dtcc+1:dtcc+4));
                   case 'baseB'
                       bsB(blcc) = mean(colH(dtcc+1:dtcc+4));
                   case 'topA'
                       tpA(blcc) = mean(colH(dtcc+1:dtcc+4));
                   case 'topB'
                       tpB(blcc) = mean(colH(dtcc+1:dtcc+4));
                   % Case Septum/Side read 11 excel rows of data
                   case 'sepA'
                       spA(blcc) = mean(colH(dtcc+1:dtcc+11));
                   case 'sepB'
                       spB(blcc) = mean(colH(dtcc+1:dtcc+11));
                  case 'sideA'
                       sdA(blcc) = mean(colH(dtcc+1:dtcc+11));
                  case 'sideB'
                       sdB(blcc) = mean(colH(dtcc+1:dtcc+11));
               end % switch/case end statement.
           else 
               if flg1 == length(colA)    % check for error.
                 clamp_er(blcc,tbcc)= 1;  % set error label.
               end
            end % terminates if statement line 43.
       end % terminates for loop dtcc
    end % terminates for loop blcc
end % terminates for loop tbcc
% ************************************************************************
% Section (2): Edit data files prior to computing deviations due to 
% additional data from coil remeasurement.
% Clamp locations 38-48: coil was shimmed at 42-45 and remeasured 38-48.
% Clamp locations 75-82: coil was shimmed ONLY on side A at clamp 78 but
% remeasured 75-82.
% ************************************************************************
% Load remeasurment data for clamps 38-48 from Excel file.
% Rows of remeasurement_delta are r1=topA, r2=topB, r3=sideA, r4=sideB.
remeasurment_deltas1=zeros(4,11); % Initializing.
remeasurment_deltas1=xlsread(strcat(s2,s1),'Measurement Edits','A62:K65');
% *** NOTE: Data in the existing Matlab files has been shifted +1 in the sequence.
% *** Thus, Excel remeasurement data clamp 38 is sequence 39, 40 is 41 etc.
% Apply changes for clamps 38-48 (sequence 39-49).
for ec1=39:49
% Adjust baseA/B data clamps 43-46 (sequence 44-47) by shim thickness 0.04 inch.
     if ec1>=43 & ec1<=46
        bsA(ec1)=bsA(ec1)+0.04;
        bsB(ec1)=bsB(ec1)+0.04;
     end
% Adjust topA/B data clamps 38-48 (sequence 39-49) by the remeasurment delta.     
     tpA(ec1)=tpA(ec1)+remeasurment_deltas1(1,ec1-38);
     tpB(ec1)=tpB(ec1)+remeasurment_deltas1(2,ec1-38);
% Adjust sideA/B data clamps 38-48 (sequence 39-49) by the remeasurment delta. 
    sdA(ec1)=sdA(ec1)+remeasurment_deltas1(3,ec1-38);
    sdB(ec1)=sdB(ec1)+remeasurment_deltas1(4,ec1-38); 
end
%-------------------------------------------------------------------------
% *** NOTE: Data in the Matlab files has been shifted +1 in the sequence.
% *** So for example, Excel remeasurement data clamp 78 is sequence 79.
% Apply changes for clamps 75-82 (sequence 76-83).
for ec1=76:83
% Adjust baseA side data clamp 78 (sequence 79) by shim thickness 0.04 inch.
     if ec1==79
        bsA(ec1)=bsA(ec1)+0.04;
        bsB(ec1)=bsB(ec1)+0.04;
      end
end

% ************************************************************************
% Section (3): Calculate Current Center deviations.
% Deviations are calculated locally (lateral (dx) & Radial (dy)) at each
% clamp location.
% ************************************************************************
% Calculate Current Center deviations locally (lateral (dx) & Radial (dy)).
dy = (bsA+bsB+tpA+tpB)/4;
% NOTE: for lateral (dx) value, septum A & side A values require negation before
% summing in order to properly account for sign/axis direction.
dx = ((sdA*-1)+(spA*-1)+sdB+spB)/4;
% ************************************************************************
% Section (3): Compare Results with Art Brook's Results.
% ************************************************************************
validation_summary=zeros(length(dx),8);
arts = xlsread(strcat(s2,s1),'Arts Ic deltas','B2:C97');
comp(:,1)=abs(dx-arts(:,1));
comp(:,2)=abs(dy-arts(:,2));
smsqcmp=sqrt(comp(:,1).^2+comp(:,2).^2);
%comp(:,1)=((round(abs(dx-arts(:,1))*1e8))/1e8);
%comp(:,2)=((round(abs(dy-arts(:,2))*1e8))/1e8);
skip_cells=sum(clamp_er,2);
for vcc=1:length(dx)
    clampn(vcc)=vcc;
    if skip_cells(vcc)==0
        validation_summary(vcc,:) = [clampn(vcc) dx(vcc) arts(vcc,1) dy(vcc) arts(vcc,2) comp(vcc,1) comp(vcc,2) smsqcmp(vcc)];
    end 
end
figure(1)
plot(validation_summary(:,8))


max_method_deltas=sortrows(validation_summary,-8);
maxdx=max(validation_summary(:,6));
maxdy=max(validation_summary(:,7));
maxSMSQ=max(validation_summary(:,8));

% ************************************************************************
% Reset Directory to Initial Matlab Directory
cd(working_matlab_dir);

% Display algorithm run time.
msgbox(strcat(num2str(toc/60),' minutes'),'Duration')