% FIELDMAPCOMP.M - Compare CST and APL field maps for RFQ
%

    clear all ; close all ;
    cd(fullfile(getgptdir,'fieldmapcomp')) ;

%    cstfile = ('/Volumes/FATSWAP/matching.txt') ;
    cstfile = fullfile(getgptdir,'fieldmapcomp','matching.txt') ;
    aplfile = fullfile(getgptdir,'rfq_apl','rfq_full','rfqmesh-full.txt') ;
    aplmap = [] ;
    cstmap = [] ;
    fieldmap = [] ;
    savefigs = 1 ;

%    Open Alan's field map for reading and read in first line
%    Columns are in the following order: z y x Ex Ey Ez Bx By Bz 

    fid = fopen(aplfile) ;
    disp(['APL fid = ' num2str(fid)]) ;

%    headdat = textscan(fid, '%s',1,'delimiter','EOL') ;
    headdat = textscan(fid, '%s %s %s %s %s %s %s %s %s',1) ;

%    Read in data, close file and add to APL field map array

    apldat = textscan(fid, '%n %n %n %n %n %n %n %n %n') ;
    fclose(fid) ;

    aplmap(:,1) = apldat{3}.*1e3 ;
    aplmap(:,2) = apldat{2}.*1e3 ;
    aplmap(:,3) = apldat{1}.*1e3 ;
    aplmap(:,4) = apldat{4} ;
    aplmap(:,5) = apldat{5} ;
    aplmap(:,6) = apldat{6} ;

%    Select required field map points

    [aplrfqinds] = find(aplmap(:,3) > 21.8) ;
    aplmap(aplrfqinds,:) = [] ;
    [aplrfqinds] = find((aplmap(:,3) == 0) | (aplmap(:,3) == 5.45) | (aplmap(:,3) == 10.9) | ...
        (aplmap(:,3) == 16.35) | (aplmap(:,3) == 21.8)) ;
    fieldmap = aplmap(aplrfqinds,:) ;

%    Open CST field map for reading and read in first line
%    Columns are in the following order: x y z Ex Ey Ez

    fid = fopen(cstfile) ;
    disp(['CST fid = ' num2str(fid)]) ;

    headdat = textscan(fid, '%s',5,'delimiter','EOL') ;

%    Read in data, close file and add to CST field map array

    cstdat = textscan(fid, '%n %n %n %n %n %n') ;
    fclose(fid) ;

    for i = 1:6
        cstmap(:,i) = cstdat{i} ;
    end

%    Select required field map points

    [cstrfqinds] = find((cstmap(:,1) < -3.5) | (cstmap(:,1) > 3.5) | (cstmap(:,2) < -3.5) | (cstmap(:,2) > 3.5)) ;
    cstmap(cstrfqinds,:) = [] ;

    yvals = ( round( [-3.5:0.7:3.5].*10 ) )/10 ;
    zvals = ( round( [0:5.45:21.8].*100 ) )/100 ;
    startind = 1 ;

    for i = 1:length(zvals)
        for j = 1:length(yvals)
            cstinds = find((cstmap(:,2) == yvals(j)) & (cstmap(:,3) == zvals(i))) ;
            fieldmap(startind:startind+10,:,2) = cstmap(cstinds,:) ;
            startind = startind + 11 ;
        end
    end

%    Plot data

    fieldmap(:,7,:) = sqrt( fieldmap(:,4,:).^2 + fieldmap(:,5,:).^2 ) ;
    figure, hold on ;
    plot(fieldmap(:,7,1)./1e6,'-b') ; plot(fieldmap(:,7,2)./1e6,'-r') ;
    xlabel(['Field data point']), ylabel(['Transverse {\it E}-field magnitude (MV/m)']), ...
        title(['\bf Transverse {\it E}-field magnitude for all data points']) ;
    axis tight ;
    hold off ; jpgplot ;
    if savefigs, figsave(gcf,'transfieldcomp_allpoints.jpg','a4') ; end ;

    fieldmap(:,8,:) = sqrt( fieldmap(:,4,:).^2 + fieldmap(:,5,:).^2 + fieldmap(:,6,:).^2 ) ;
    figure, hold on ;
    plot(fieldmap(:,8,1)./1e6,'-b') ; plot(fieldmap(:,8,2)./1e6,'-r') ;
    xlabel(['Field data point']), ylabel(['Full {\it E}-field magnitude (MV/m)']), ...
        title(['\bf Full {\it E}-field magnitude for all data points']) ;
    axis tight ;
    hold off ; jpgplot ;
    if savefigs, figsave(gcf,'fullfieldcomp_allpoints.jpg','a4') ; end ;

    for i = 1:5
        figure, hold on ;
        startind = ((i-1)*121) + 1 ;
        endind = i*121 ;
        if max(fieldmap(startind:endind,7,1)) == 0 || max(fieldmap(startind:endind,7,2)) == 0
            scalefac = 1 ;
        else
            scalefac = ( max(fieldmap(startind:endind,7,2)) )./( max(fieldmap(startind:endind,7,1)) ) ;
        end
        quiver(fieldmap(startind:endind,1,1),fieldmap(startind:endind,2,1), ...
            fieldmap(startind:endind,4,1),fieldmap(startind:endind,5,1),'-b','filled') ;
        quiver(fieldmap(startind:endind,1,2),fieldmap(startind:endind,2,2), ...
            -fieldmap(startind:endind,4,2),-fieldmap(startind:endind,5,2),scalefac,'-r','filled') ;
        xlabel(['{\it x}-position (mm)']), ylabel(['{\it y}-position (mm)']), ...
            title(['\bf Transverse {\it E}-field at {\it z}=' num2str(zvals(i)) ' mm']) ;
        axis([-5 5 -5 5]) ;
        axis equal ; axis square ;
        hold off ; jpgplot ;
        if savefigs, figsave(gcf,['Efield_quivplot_' num2str(zvals(i)) 'mm.jpg'],'screen') ; end ;
    end

    return

    