function out = return_profile_plots(in,consts)

% This code generates GMT plots showing production profiles
% from Cl, Ca, K targets. 
% Needs to have scaling factor part removed from this and put upstream

if nargin < 2
    load Cl36_consts;
end


% GMT initialize

formatstr = '--FONT_TITLE=12p,3,black --MAP_ANNOT_OFFSET_PRIMARY=0.1i --MAP_TICK_LENGTH_PRIMARY=-0.05i --MAP_FRAME_PEN=0.5p --FONT_ANNOT_PRIMARY=10p,2,black --FONT_LABEL=12p,3,black';

ymax = 100;
ysp = 10;

out = {};

% Repeat for all samples in results...

for a = 1:in.numcls
    % Get filename
    uniqueid = int2str(round(rand(1)*1e6));
    fnamestub = ['profile36_v3_' uniqueid];
    
    % Get data
    
    p = in.profile{a};
    plot_z = p.zz;
    
    % Profile is actually in normalized units for spallation. Compute
    % actual spallation using St scaling
    plot_Psp_Ca = p.P36sp_Ca.*consts.Psp_0_Ca_St.*in.n.sf_St(a);
    plot_Psp_K = p.P36sp_K.*consts.Psp_0_K_St.*in.n.sf_St(a);
    plot_Psp_Ti = p.P36sp_Ti.*consts.Psp_0_Ti_St.*in.n.sf_St(a);
    plot_Psp_Fe = p.P36sp_Fe.*consts.Psp_0_Fe_St.*in.n.sf_St(a);
    
    % For anything to do with spallogenic thermal/epithermal neutron 
    % capture, apply St scaling factor as well
    plot_Psp_n = consts.P_f_0_St.*(p.P36sp_th_Cl + p.P36sp_eth_Cl).*in.n.sf_St(a);
    
    % Muon-related stuff already calculated
    plot_Pmu_n = p.P36mu_eth_Cl + p.P36mu_th_Cl;
    plot_Pmu_K = p.P36_neg_K + p.P36_fast_K;
    plot_Pmu_Ca = p.P36_neg_Ca + p.P36_fast_Ca;
    
    % Plotting by target
    plot_PCa = plot_Psp_Ca + plot_Pmu_Ca;
    plot_PK = plot_Psp_K + plot_Pmu_K;
    plot_PTi = plot_Psp_Ti;
    plot_PFe = plot_Psp_Fe;
    plot_PCl = plot_Psp_n + plot_Pmu_n;
    
    % Total production
    plot_P = plot_Psp_Ca + plot_Psp_K + plot_Psp_Fe + plot_Psp_Ti + plot_Psp_n + plot_Pmu_n + plot_Pmu_K + plot_Pmu_Ca;
        
    % ------------ Determine x axis limit -------------------
    
    maxx = max(plot_P);

    possx = [10 20 50 100 200 500 1000 2000];
    possxi = [2 5 10 20 50 100 200 500];

    uselim = min(find(possx > maxx));
    if isempty(uselim); uselim = length(possx); end

    xmax = possx(uselim);
    xsp = possxi(uselim);
    
    % ------------ Initialize GMT ------------------
    
    fname = [consts.scratchdir fnamestub];
    gmtname = [fname '.gmt'];
    psname = [fname '.ps'];
    pngname = [fname '.png'];
    
    % open file for write
    fid = fopen(gmtname,'w');

    % Create GMT script
    % Define format strings

    lxs = ['P36 (atoms/g/yr)'];
    lys = ['Depth (g/cm2)'];


    rstring = ['-R0/' sprintf('%0.0f',xmax) '/0/' sprintf('%0.0f',ymax)];
    jstring = '-JX3i/-3i';
    bstring = ['-BWeSn+t"' in.s.sample_name{a} '" -Bpxa' sprintf('%0.0f',xsp) 'g' sprintf('%0.0f',xsp) '+l"' lxs '" -Bpya' sprintf('%0.0f',ysp) 'g' sprintf('%0.0f',ysp) '+l"' lys '"'];
    
    % Tell shell path to GMT if necessary, this is for GB's Mac
    if consts.isLocal == 1
        fprintf(fid,'%s\n\n','PATH=$PATH:/sw/bin');
        fprintf(fid,'%s\n\n','export PATH');
    end
    
    % GMT basemap
    % Make the GMT basemap
    fprintf(fid,'%s\n','# ----------- Draw the plot axes ------------');

    temp_string = ['gmt psbasemap ' bstring ' ' jstring ' ' rstring ' -K -P ' formatstr ' > ' psname];
    fprintf(fid,'%s\n',temp_string);
    
    % Plot profiles
    % Color coding: black, total production
    % Color by target: green, Cl; dark blue, Ca, light blue, K, orange, Ti,
    % red, Fe
    
    fprintf(fid,'%s\n','# ----------- Plot Cl production ------------');
    
    temp_string = ['gmt psxy ' rstring ' ' jstring ' -P -O -K -W1p,0/175/0 << EOD >> ' psname];
    fprintf(fid,'%s\n',temp_string);
    fprintf(fid,'%0.2e %0.2e\n',[plot_PCl' plot_z']');
    fprintf(fid,'%s\n','EOD');
    
    fprintf(fid,'%s\n','# ----------- Plot Ca production ------------');
    
    temp_string = ['gmt psxy ' rstring ' ' jstring ' -P -O -K -W1p,0/0/175 << EOD >> ' psname];
    fprintf(fid,'%s\n',temp_string);
    fprintf(fid,'%0.2e %0.2e\n',[plot_PCa' plot_z']');
    fprintf(fid,'%s\n','EOD');
    
    fprintf(fid,'%s\n','# ----------- Plot K production ------------');
    
    temp_string = ['gmt psxy ' rstring ' ' jstring ' -P -O -K -W1p,150/150/255 << EOD >> ' psname];
    fprintf(fid,'%s\n',temp_string);
    fprintf(fid,'%0.2e %0.2e\n',[plot_PK' plot_z']');
    fprintf(fid,'%s\n','EOD');
    
    fprintf(fid,'%s\n','# ----------- Plot Ti production ------------');
    
    temp_string = ['gmt psxy ' rstring ' ' jstring ' -P -O -K -W1p,255/165/0 << EOD >> ' psname];
    fprintf(fid,'%s\n',temp_string);
    fprintf(fid,'%0.2e %0.2e\n',[plot_PTi' plot_z']');
    fprintf(fid,'%s\n','EOD');
    
    fprintf(fid,'%s\n','# ----------- Plot Fe production ------------');
    
    temp_string = ['gmt psxy ' rstring ' ' jstring ' -P -O -K -W1p,150/0/0 << EOD >> ' psname];
    fprintf(fid,'%s\n',temp_string);
    fprintf(fid,'%0.2e %0.2e\n',[plot_PFe' plot_z']');
    fprintf(fid,'%s\n','EOD');
    
    fprintf(fid,'%s\n','# ----------- Plot total production ------------');
    
    temp_string = ['gmt psxy ' rstring ' ' jstring ' -P -O -K -W1p,0 << EOD >> ' psname];
    fprintf(fid,'%s\n',temp_string);
    fprintf(fid,'%0.2e %0.2e\n',[plot_P' plot_z']');
    fprintf(fid,'%s\n','EOD');
    
    fprintf(fid,'%s\n','# ------------ Plot legend ------------------');
    fprintf(fid,'%s\n',['gmt pslegend -O ' jstring ' -Dx2.15i/0.1i+w0.75i/1.15i -F+gwhite << EOD >> ' psname]);
    fprintf(fid,'%s\n','S 0.1i - 0.25i - 1p,0 0.3i Total');
    fprintf(fid,'%s\n','S 0.1i - 0.25i - 1p,0/175/0 0.3i Cl');
    fprintf(fid,'%s\n','S 0.1i - 0.25i - 1p,150/150/255 0.3i K');
    fprintf(fid,'%s\n','S 0.1i - 0.25i - 1p,0/0/175 0.3i Ca');
    fprintf(fid,'%s\n','S 0.1i - 0.25i - 1p,255/165/0 0.3i Ti');
    fprintf(fid,'%s\n','S 0.1i - 0.25i - 1p,150/0/0 0.3i Fe');
    fprintf(fid,'%s\n','EOD');
    
    
    % close file 
    fclose(fid);
    
    if consts.isLocal == 0
        % Then you are running on a linux system...
        % Execute the script you just made
        system(['chmod a+x ' gmtname]);
        system(gmtname);

        % convert the ps file to a png using ImageMagick...
        system(['convert ' psname ' -trim -transparent white ' pngname]);       
    else
        % On GB's Mac
        disp(['sh ' gmtname]);
        disp(['open ' psname]);
        % Running the convert command doesn't seem to work. 
    end
    
    % record plot name
    out{end+1} = fnamestub;
    
end % close main loop

% return list of plotnames
    