Commit fe219ed9 authored by Lore, Jeremy's avatar Lore, Jeremy
Browse files

LP ELM filtering added

parent b1a4b46f

matlab/bfield_library_jdl/D3D/plot_ir_profile.m

+4 −4
Original line number Diff line number Diff line
@@ -21,13 +21,13 @@ elmWin = [0.8,0.95]; 
sparseifyLines = 20;    % Step for time trace indices, only in plotting

crangeElmPlot = 2;  % Colorbar control for ELM cycle plot: 1 = [0,1], else tight to elmWin set above



elmThreshNorm = 0.3;   % In normalized magnitude for ELM filter

elmThreshNorm = 0.2;   % In normalized magnitude for ELM filter

dtThreshMS = 1;  % Discard "ELM" spikes repeated in this threshold (ms)

% Select ir radius to find ELMS (in dR from Outer SP)

dROutFindELMCM = -10;



WRITE_FILE = 0;

MAKE2DPLOTS = 1;

MAKE2DPLOTS = 0;



dROuterWantRangeCM = [-10,20];

dRInnerWantRangeCM = [-5,8];
@@ -42,7 +42,7 @@ nInnerBinIR = 12; % Number of radial bins in dR range 


%% Load file(s)

IR = get_IR_data(fname,tWinMS);

IR = filterDataElms(IR,dROutFindELMCM,elm_file,elmThreshNorm,dtThreshMS);

IR = tagDataElms(IR,dROutFindELMCM,elm_file,elmThreshNorm,dtThreshMS);

if MAKE2DPLOTS

    plotIR2D(IR,dROuterWantRangeCM,dRInnerWantRangeCM);

end
@@ -196,7 +196,7 @@ end 






function IR = filterDataElms(IR,dROutFindELMCM,elm_file,elmThreshNorm,dtThreshMS)

function IR = tagDataElms(IR,dROutFindELMCM,elm_file,elmThreshNorm,dtThreshMS)



%% ELM Filter

if ~isempty(elm_file)

matlab/bfield_library_jdl/D3D/process_LP.m

+100 −43
Original line number Diff line number Diff line
@@ -3,9 +3,15 @@ clearvars; 
% Channels 1:8 are on shelf, increasing in radius

% 9:19 are inner div, moving down center stack then out in radius



tWinMS = [3100,4000];

gfile_name = 'C:\Users\jjl\Dropbox (ORNL)\DIII-D\Qprl experiment\power13mw\g174310.03500_153';

tWinMS = [3200,4000];

elmWin = [0.3,0.7];

crangeElmPlot = 2;  % Colorbar control for ELM cycle plot: 1 = [0,1], else tight to elmWin set above

elmThreshNorm = 0.1;

elmTimeNoRepeatMS = 1;

signalNamePlot = 'heatflux';  % Just used to make figures

% gfile_name = 'C:\Users\jjl\Dropbox (ORNL)\DIII-D\Qprl experiment\power13mw\g174310.03500_153';

fileNameLP = 'C:\Users\jjl\Dropbox (ORNL)\DIII-D\Qprl experiment\1743XX\LP_174306.mat';

elmFile = 'C:\Users\jjl\Dropbox (ORNL)\DIII-D\Qprl experiment\1743XX\fs_174306_2000_5000.mat';



[filepath,name,ext] = fileparts(fileNameLP);

outfile = fullfile(filepath,[name,'_processed',ext]);
@@ -15,7 +21,7 @@ chanShelf = 1:8; 
chanInner = 9:19;



LP = load(fileNameLP); LP = LP.lp;

g = readg_g3d(gfile_name);

% g = readg_g3d(gfile_name);





%% Filter on tWin
@@ -23,30 +29,95 @@ LP = getLPDataTWin(LP,tWinMS); 
LP = divideLPShelfInner(LP,chanShelf,chanInner);

LP = makeLPRadialProf(LP);



elmData = load(elmFile); tElm = double(elmData.fs.t00); sElm = double(elmData.fs.fs00);

iStartElmFilt = find(tElm > tWinMS(1),1,'first');

iEndElmFilt = find(tElm > tWinMS(2),1,'first');

Spikes = find_spikes(tElm(iStartElmFilt:iEndElmFilt),sElm(iStartElmFilt:iEndElmFilt),elmThreshNorm,elmTimeNoRepeatMS,1);



LP = filterLPElms(LP,Spikes,elmWin);



figure; hold on; box on; grid on; set(gcf,'color','w');

for i = 1:length(LP.Shelf.chan)

    thisChan = LP.Shelf.chan{i};

%     plot(thisChan.delrsepout,thisChan.temp,'.')    

    plot(thisChan.dLSepOut,thisChan.jsat,'x')    

    plot(thisChan.dLSepOut,thisChan.(signalNamePlot),'x')    

    xlabel('dL_{sep} Outer')

    ylabel('j_{sat} ()')

    ylabel(signalNamePlot)

    title('Shelf')

end

plot(LP.Profiles.Shelf.dLSepOut,medflt1d_jl(LP.Profiles.Shelf.(signalNamePlot),100))

% plot(LP.ProfilesELMFiltered.Shelf.dLSepOut,LP.ProfilesELMFiltered.Shelf.(signalNamePlot),'ko')

plot(LP.ProfilesELMFiltered.Shelf.dLSepOut,medflt1d_jl(LP.ProfilesELMFiltered.Shelf.(signalNamePlot),100),'k')













num_col = 1024;



if crangeElmPlot == 1

    crange = [0,1];

else

    crange = elmWin;

end

icuse = interp1(linspace(crange(1),crange(2),num_col),1:num_col,LP.Profiles.Shelf.periodELMCycle);



figure; hold on; box on; grid on; set(gcf,'color','w');

iBad = LP.Profiles.Shelf.periodELMCycle == -1;

scatter(LP.Profiles.Shelf.dLSepOut(~iBad),LP.Profiles.Shelf.(signalNamePlot)(~iBad),12,LP.Profiles.Shelf.periodELMCycle(~iBad),'filled')

colorbar

xlabel('dL_{sep} Outer')

ylabel(signalNamePlot)

title('Shelf')

colormap(colorflipper(num_col,'plasma'));





figure; hold on; box on; grid on; set(gcf,'color','w');

for i = 1:length(LP.Inner.chan)

    thisChan = LP.Inner.chan{i};

%     plot(thisChan.delrsepin,thisChan.temp,'.')    

    plot(thisChan.dLSepIn,thisChan.jsat,'x')   

    plot(thisChan.dLSepIn,thisChan.(signalNamePlot),'x')   

    xlabel('dL_{sep} Inner')

    ylabel('j_{sat} ()')

    ylabel(signalNamePlot)

    title('Inner')

end

plot(LP.ProfilesELMFiltered.Inner.dLSepIn,medflt1d_jl(LP.ProfilesELMFiltered.Inner.(signalNamePlot),100),'k')



save(outfile,'LP');





function LP = filterLPElms(LP,Spikes,elmWin)



locations = {'Shelf','Inner'};

for i = 1:length(locations)

    thisLocationName = locations{i};



    %% Tag time in ELM cycle

    thisTimeSeries = LP.Profiles.(thisLocationName).time;

    LP.Profiles.(thisLocationName).periodELMCycle = -1*ones(size(thisTimeSeries));

    

    for j=1:length(thisTimeSeries)

        ind2 = find(Spikes.times(1:end-1) <= thisTimeSeries(j) & Spikes.times(2:end) > thisTimeSeries(j));

        if ~isempty(ind2)

            dte2 = Spikes.times(ind2+1)-Spikes.times(ind2);

            dts2 = thisTimeSeries(j) - Spikes.times(ind2);

            LP.Profiles.(thisLocationName).periodELMCycle(j) = dts2/dte2;

        end

    end    

    

    % Filter using elmWin

    fieldNames = fieldnames(LP.Profiles.(thisLocationName));

    inElmCycle = find(LP.Profiles.(thisLocationName).periodELMCycle >= elmWin(1) & LP.Profiles.(thisLocationName).periodELMCycle <= elmWin(2));

    for k = 1:length(fieldNames)

        thisFieldName = fieldNames{k};        

        LP.ProfilesELMFiltered.(thisLocationName).(thisFieldName) = LP.Profiles.(thisLocationName).(thisFieldName)(inElmCycle);

    end

    

end



end



function LP = divideLPShelfInner(LP,chanShelf,chanInner)



iCount = 1;
@@ -84,8 +155,8 @@ for j = 1:length(locations) 
            if length(thisChan.(thisFieldName)) == nD

                LP.Profiles.(thisLocation).(thisFieldName)(theseInds) = thisChan.(thisFieldName);

            end            

            indexOffset = indexOffset + nD;

        end

        indexOffset = indexOffset + nD;

    end

        

end
@@ -116,32 +187,18 @@ function LPOut = getLPDataTWin(LP,tWinMS) 
    end

    

    keepInds = find(thisChan.time >= tWinMS(1) & thisChan.time <= tWinMS(2));

        LPOut.chan{i} = thisChan;

    fieldNames = fieldnames(thisChan);

    for k = 1:length(fieldNames)

        thisFieldName = fieldNames{k};

        if length(thisChan.(thisFieldName)) == length(thisChan.time)

            LPOut.chan{i}.(thisFieldName) = double(thisChan.(thisFieldName)(keepInds));

        else

            LPOut.chan{i}.(thisFieldName) = thisChan.(thisFieldName);

        end        

    end

    LPOut.chan{i}.ntimes = length(keepInds);

        LPOut.chan{i}.time = double(thisChan.time(keepInds));

        LPOut.chan{i}.temp = double(thisChan.temp(keepInds));

        LPOut.chan{i}.temp_err = double(thisChan.temp_err(keepInds));

        LPOut.chan{i}.dens = double(thisChan.dens(keepInds));

        LPOut.chan{i}.dens_err = double(thisChan.dens_err(keepInds));

        LPOut.chan{i}.pot = double(thisChan.pot(keepInds));

        LPOut.chan{i}.pot_err = double(thisChan.pot_err(keepInds));

        LPOut.chan{i}.isat = double(thisChan.isat(keepInds));

        LPOut.chan{i}.jsat = double(thisChan.jsat(keepInds));

        LPOut.chan{i}.angle = double(thisChan.angle(keepInds));

        LPOut.chan{i}.area = double(thisChan.area(keepInds));

        LPOut.chan{i}.psin = double(thisChan.psin(keepInds));

        LPOut.chan{i}.delrsepin = double(thisChan.delrsepin(keepInds));

        LPOut.chan{i}.delrsepout = double(thisChan.delrsepout(keepInds));

        LPOut.chan{i}.delzsepin = double(thisChan.delzsepin(keepInds));

        LPOut.chan{i}.delzsepout = double(thisChan.delzsepout(keepInds));

        LPOut.chan{i}.csq = double(thisChan.csq(keepInds));

        LPOut.chan{i}.res_err = double(thisChan.res_err(keepInds));

        LPOut.chan{i}.heatflux = double(thisChan.heatflux(keepInds));

        

    LPOut.chan{i}.dLSepOut = sqrt( LPOut.chan{i}.delrsepout.^2 + LPOut.chan{i}.delzsepout.^2);

    LPOut.chan{i}.dLSepIn  = sqrt( LPOut.chan{i}.delrsepin.^2  + LPOut.chan{i}.delzsepin.^2 );

end



end
 


% plot(LP.hf_r_rsep+Rsep,medflt1d_jl(LP.hf,5),'rx','HandleVisibility','off')
 
 No newline at end of file