Add some older DIII-D expt data routines

function elm=check_elm_timing(elm_file,twin)

% clearvars;

myelm = load(elm_file);

% load(fs_file);

% twin = [3000,4900];

elm_norm_thresh = 0.4;

elm_dt_thresh = 10;

igood = find(myelm.elm.telmpeak >= twin(1) & myelm.elm.telmpeak <= twin(2));

telm = myelm.elm.telmpeak(igood);

elmpeak = myelm.elm.elmpeak(igood);

igood_fs = find(myelm.elm.t >= twin(1) & myelm.elm.t <= twin(2));

tfs = myelm.elm.t(igood_fs);

fs = myelm.elm.d(igood_fs);

emean = mean(elmpeak);

elmpeak_norm = elmpeak/max(elmpeak);

ielm_keep = find(elmpeak_norm >= elm_norm_thresh);

tkeep = telm(ielm_keep);

% dt filter

dtkeep = diff(tkeep);

ikeep_dt = find(dtkeep >= elm_dt_thresh);

ielm_keep = ielm_keep(ikeep_dt);

% Final results

tkeep = tkeep(ikeep_dt);

nkeep = length(tkeep);

fs_at_elm = interp1(tfs,fs,tkeep);

figure; hold on;

plot(tfs,fs./max(fs),'k')

plot(telm,elmpeak_norm,'r.')

% plot(telm(ielm_keep),elmpeak_norm(ielm_keep),'ro')

plot(tkeep,fs_at_elm./max(fs),'ro')

% plot(telm,emean*ones(size(telm)),'k')

ylabel('Norm. filterscope')

xlabel('Time [ms]')

title('ELM timing')

legend('fs','elm marks','selected')

elm.telm = tkeep;

elm.nelm = nkeep;

clearvars;

gfile_name = 'C:\Work\DIII-D\APS 2016\g156855.04526_694';

g = readg_g3d(gfile_name);

files{1} = 'C:\Work\DIII-D\APS 2016\profs_156855_4500_awlhm.mat'; mytitle{1} = '156855';

files{2} = 'C:\Work\DIII-D\APS 2016\profs_156859_4500_awlhm.mat'; mytitle{2} = '156859';

% files{2} = 'C:\Work\DIII-D\APS 2016\profs_156861_4500_awlhm.mat'; mytitle{2} = '156861';

files{3} = 'C:\Work\DIII-D\APS 2016\profs_156867_4500_awlhm.mat'; mytitle{3} = '156867';

files{4} = 'C:\Work\DIII-D\APS 2016\profs_156871_4500_awlhm.mat'; mytitle{4} = '156871';

for i = 1:length(files)

    myprofs{i} = load(files{i});

end

psiNfit = linspace(0.95,1.05,1000);

psiN_core_SOLPS = 0.9;

tewants = [63,63,85,83];

shifts = zeros(size(tewants));

for i = 1:length(files)

    tewant = tewants(i)/1000;

    profs = myprofs{i}.profs;

    tefit = evaluate_tanh_fit(profs.tetanh,psiNfit);    

    Te_sep_noshift = evaluate_tanh_fit(profs.tetanh,1);

    shifts(i) = 1-interp1(tefit,psiNfit,tewant);

    Te_check = evaluate_tanh_fit(profs.tetanh,1-shifts(i));

    [rs,zs] = calc_RZ_at_psiN_theta(g,[1,1-shifts(i)],0);

    real_shift = -diff(rs);

    % Find n at core SOLPS boundary

    nfit = evaluate_tanh_fit(profs.ntanh,psiN_core_SOLPS-shifts(i));

    fprintf('\n')

    fprintf('%s\n',mytitle{i})

    fprintf('Te_sep, no shift = %f [eV]\n',Te_sep_noshift*1000)

    fprintf('Want Te_sep      = %f [eV]\n',tewant*1000)

    fprintf('Got Te_sep       = %f [eV]\n',Te_check*1000)

    fprintf('psiN shift       = %f [eV]\n',shifts(i)*1000)

    fprintf('dR shift [mm]    = %f [eV]\n',real_shift*1000)

    fprintf('ne at core bdry  = %f [e20] \n',nfit)

%     fprintf('sqrt(pshift)[~mm]= %f [eV]\n',sqrt(shifts(i)*1000))

end

 No newline at end of file

function spikes = find_spikes(t,d,thresh,dt_thresh,MEDFILT)

if MEDFILT

d=d-medflt1d_jl(d,floor(length(t)*0.1),'zero');

end

med = median(d);

sig_norm = (d - med)./max(d-med);

% Find peaks

right_neighbor = [0;sig_norm(1:end-1)];

left_neighbor = [sig_norm(2:end);0];

p_ind = find((sig_norm > right_neighbor) & (sig_norm > left_neighbor) & (sig_norm > thresh));

% eliminate multiples

dt_ind = diff(t(p_ind));

max_loop = 10;

iloop = 0;

while any(dt_ind < dt_thresh)

    i = 1;

    while i < length(dt_ind) + 1

        if dt_ind(i) < dt_thresh %&& dt_ind(i+1) < dt_thresh

            if sig_norm(p_ind(i)) > sig_norm(p_ind(i+1))

                p_ind(i+1) = [];

                i = i - 1;

            else

                p_ind(i) = [];

                i = i - 1;

            end

        end

        dt_ind = diff(t(p_ind));

        i = i+1;

    end

    iloop = iloop + 1;

    if iloop > max_loop

        disp(dt_ind)

        error('too many iterations')

    end    

end

figure; hold on; box on;

subplot(2,1,1); hold on; box on;

plot(t,d)

yline(med,'k')

subplot(2,1,2); hold on; box on;

ylabel('Heat flux [MW/m^2]')

xlabel('Time [ms]')

plot(t,sig_norm)

plot(t(p_ind),sig_norm(p_ind),'ro')

yline(thresh,'k')

ylabel('Normalized signal')

xlabel('Time [ms]')

title('ELM timing (IR)')

spikes.time_inds = p_ind;

spikes.times = t(p_ind);

spikes.val = d(p_ind);

spikes.val_norm = sig_norm(p_ind);

 No newline at end of file

clearvars;

fname = 'C:\Work\DIII-D\APS 2016\irtv01_156855_2000_4900.mat'; shot = 156855; elm_file = 'elm_156855_fs6middaf.mat';

load(fname);

thresh = 0.5;

dt_thresh = 10; %ms

twin = [4200,4800];

% twin = [3200,3800];

% elm_win = [0.8,0.95];

elm_win = [0.1,0.9];

DIV_PLOT = 1; % 1 = outer

if DIV_PLOT == 1

    drwant_min = -10;

    drwant_max = 20;

else

    error('set up inner divertor')

end

% RAW

time = double(ir.time);

data = double(ir.data);

drout = double(ir.drout);

drin = double(ir.drin);

% TIME SELECT

t_good = find(time >= twin(1) & time <= twin(2));

% eskeep = es(t_good);

time = time(t_good);

data = data(:,t_good);

drout = drout(:,t_good);

drin = drin(:,t_good);

drfind_elm = -10;

for i = 1:length(time)

    [dx,ix] = min(abs(drout(:,i) - drfind_elm));

    dat_1d(i) = data(ix,i);

end

sig = dat_1d;

spikes = find_spikes(time,sig,thresh,dt_thresh);

clearvars;

run_path = 'C:\Users\jjl\Dropbox (ORNL)\DIII-D\Qprl experiment\1743XX';

outFileName = 'downstream_data_174306.mat';

% From plot_ir_profile

IRfileName = 'IR_profile_174306_twin_3100_4000_elm_80_95.mat';

IR = load(fullfile(run_path,IRfileName)); IR = IR.IR;

DownStream.HeatFlux.IR = IR;

save(fullfile(run_path,outFileName),'DownStream')

 No newline at end of file