Plot with 2 corresponding x-axes

This question may seem like it has been answered in the forums but I believe there is a subtlety that has not been directly answered. I want to plot the optical absorption coefficient of a material as a function of both photon wavelength (bottom x-axis) and energy (top x-axis) similar to this figure . There are two problems with this task that I have not found sufficiently addressed. The first is to ensure that all of the points in the bottom axis line up with the points in the top axis (we need to some how link the two x-axes). Second, because of how energy and wavelength are related, one axis will be ascending while the other will be descending (Matlab doesn't like descending x-axes). How can you make the referenced plot? Can someone recreate the reference plot using simplified data?

x1 = 1:10;
x2 = 1./x1;
y = -x1 + 10;

Answer 1

[dpb]

Basics aren't difficult; I'll have to consider precisely how to do the upper axis ticks; I've got another appointment in town shortly so no time to piddle but I think you'll have to use linear scale and then compute the percentage for each major/minor tick because it's nonlinear and there's not builtin feature...but to get the general outline isn't too bad...

figure
hAx(1)=axes;                                   % first axes
hAx(2)=axes('Position',hAx(1).Position, ...    % 2nd on top of first
            'color','none', ...                % so will show through
            'XAxisLocation','top', ...         % move XAxis to top
            'XDir','reverse');                 % and flip LR
hAx(1).XLim=[0.2 1.8];                         % first axes limits
eV=hb./hAx(1).XLim*6.242E18*1E6;               % limits of 2nd in eV from first
hAx(2).XLim=sort(eV);                          % match the limits
set(hAx,{'YScale'},{'log'},{'YLim'},{[1E3 1E8]});  % log y axes

... I don't have time to work out the algebra just now, sorry...if you don't get it to work and somebody else doesn't come along to finish it up I'll try to look in again tomorrow.

ADDENDUM

OK, carrying on from above,

fneV=@(lam) hb./lam*6.242E18*1E6;        % function for wavelength-->eV
fnWL=@(eV) hb./eV*6.242E18*1E6;          % and vice versa...convenience
hAx(2).XLim=hAx(1).XLim;                 % use same limits to match up to
lamAx=fnWL([0.7:0.1:0.9 1:5]);           % wavelengths for the wanted eV values
% draw ticks at matching wavelength, label with eV value
set(hAx(2),'XTick',sort(lamAx),'XTickLabel',fneV(lamAx))
hAx(2).FontSize=8;                       % need a little more room
hAx(1).YTick=[];                         % remove superfluous axis ticks
box on                                   % add the box around the figure

This yields the basic axes in which to draw as function of wavelength into hAx(1) that looks like

I believe the sample figure eV ticks don't quite line up at where they're supposed to be; looks to me like maybe they were approximated with a log axis instead??? Unless I made an error in one of the constants in defining the transformation equations but I don't believe I did...

This has only one really ugly artifact that is difficult to overcome in HG2; the "doubled-up" ticks on the lower x-axis; those are owing to box,'on' for hAx(2); doing that turns on the tick marks on the opposite axis from the primary one to match and since in this case the two scales don't align, the the ticks compete with each other. The only way I've imagined to solve this is to turn the box off for both axes and add yet a third to outline the plot, or draw a physical line to close the axes box. Anything done is an ugly kludge, however...

[Christopher Saltonstall]

Did you ever find a solution to getting ride of the extra tick marks?

[dpb]

The only solution is one of those outlined above...turn the box off so opposite axes ticks don't show is only way to make them disappear. Sad but true; TMW hasn't changed anything in the internals to let display the box but not draw the ticks on the opposite side.

While I've not actually done it, shouldn't be too hard to just draw a rectangle on top to simulate...

ADDENDUM:

hAx=axes;
xl=xlim; yl=ylim;
hR=rectangle('Position',[xl(1) yl(1) diff(xl) diff(yl)],'EdgeColor','r');

[Christopher Saltonstall]

what is hb in your code?

[Christopher Saltonstall]

I am not getting the same result as you even though I directly copied your code. The following is my code. Note that I have included data from the reference plot. So the results can be compared to it.

clear
close all

data = [0.35737704918032787, 92649850.48039015
    0.3819672131147541, 72475211.53588514
    0.41147540983606556, 55967196.301264346
    0.4229508196721312, 51126728.207795605
    0.4557377049180328, 47937263.23846833
    0.5016393442622952, 43803474.996365294
    0.539344262295082, 40544606.59504045
    0.5688524590163934, 34277014.89316264
    0.5950819672131147, 25461360.619704828
    0.6229508196721312, 16618285.037019765
    0.6655737704918033, 10037712.197086804
    0.7213114754098361, 6724703.556947659
    0.7950819672131147, 4169337.0144349397
    0.8836065573770491, 2722604.5977291227
    0.9557377049180328, 1997130.2156342946
    1.0344262295081967, 1523018.3932672704
    1.1049180327868853, 1255113.4937515096
    1.1852459016393442, 1034416.372625131
    1.2737704918032788, 886319.8712540427
    1.359016393442623, 779309.6059234422
    1.4540983606557378, 609961.3425615291
    1.4918032786885247, 502533.5248888627
    1.5098360655737708, 398203.214420678
    1.5196721311475412, 299601.3683763488
    1.5213114754098362, 195504.2893801002
    1.5262295081967214, 119588.7737216253
    1.5360655737704918, 65114.84416472715
    1.5475409836065577, 38316.26954492593
    1.5508196721311474, 25329.082065449114
    1.559016393442623, 18331.533361062826
    1.5688524590163935, 14713.897024069562
    1.5934426229508198, 11361.99263934681
    1.6196721311475408, 8773.786275376307
    1.6508196721311474, 6433.761487810411
    1.6754098360655738, 5032.800618487584
    1.7049180327868854, 3787.196421050415
    1.7327868852459019, 2886.942647505133
    1.7557377049180327, 2287.679320744783
    1.777049180327869, 1933.9100405245356];


%planks constant
h = 4.135e-15; %eV s
%speed of light
c = 3e8; %m/s

%absorption coef
alpha = data(:,2); %1/m

%wavelength
lambda = data(:,1); %microns

%energy
E = h*c./(lambda*1e-6);


figure(1)
plot(lambda,alpha)
xlabel('Wavelength (\mum)')
ylabel('\alpha (m^{-1})')
set(gca,'Yscale','log')

figure(2)
plot(E,alpha)
xlabel('Energy (eV)')
ylabel('\alpha (m^{-1})')
set(gca,'Yscale','log')


hb = 1;

figure(3)
%plot(lambda,alpha)
hAx(1)=axes;                                   % first axes
hAx(2)=axes('Position',hAx(1).Position, ...    % 2nd on top of first
    'color','none', ...                % so will show through
    'XAxisLocation','top', ...         % move XAxis to top
    'XDir','reverse');                 % and flip LR
hAx(1).XLim=[0.2 1.8];                         % first axes limits
eV=hb./hAx(1).XLim*6.242E18*1E6;               % limits of 2nd in eV from first
hAx(2).XLim=sort(eV);                          % match the limits
set(hAx,{'YScale'},{'log'},{'YLim'},{[1E3 1E8]});  % log y axes
fneV=@(lam) hb./lam*6.242E18*1E6;        % function for wavelength-->eV
fnWL=@(eV) hb./eV*6.242E18*1E6;          % and vice versa...convenience
hAx(2).XLim=hAx(1).XLim;                 % use same limits to match up to
lamAx=fnWL([0.7:0.1:0.9 1:5]);           % wavelengths for the wanted eV values
% draw ticks at matching wavelength, label with eV value
set(hAx(2),'XTick',sort(lamAx),'XTickLabel',fneV(lamAx))
hAx(2).FontSize=8;                       % need a little more room
hAx(1).YTick=[];                         % remove superfluous axis ticks
box on                                   % add the box around the figure

[dpb]

Hmmm....been so long! I'm pretty sure hB was an axes handle from command window that probably just typed ha=..., hb=... when first messing around. Then missed getting the references fixed when cleaning up code. Hence, I'm guessing that should have been

eV=hAx(2).XLim./hAx(1).XLim*6.242E18*1E6; 

which would be the two limits of the second axes in eV units instead of wavelength. Remember, that's from quick read after almost two years! :)

But makes sense -- need to scale the two axes to match the endpoints in the units of each to be able to place corresponding data at the right locations of each.

I have very vague recollections of this, but I didn't do any of it in the editor, only at the command line so there will be no record or script files to recover, unfortunately, at this point.

[Christopher Saltonstall]

I figured out what hb was. However, when I try to add a plot to the code, everything falls apart.

clear
close all

data = [0.35737704918032787, 92649850.48039015
    0.3819672131147541, 72475211.53588514
    0.41147540983606556, 55967196.301264346
    0.4229508196721312, 51126728.207795605
    0.4557377049180328, 47937263.23846833
    0.5016393442622952, 43803474.996365294
    0.539344262295082, 40544606.59504045
    0.5688524590163934, 34277014.89316264
    0.5950819672131147, 25461360.619704828
    0.6229508196721312, 16618285.037019765
    0.6655737704918033, 10037712.197086804
    0.7213114754098361, 6724703.556947659
    0.7950819672131147, 4169337.0144349397
    0.8836065573770491, 2722604.5977291227
    0.9557377049180328, 1997130.2156342946
    1.0344262295081967, 1523018.3932672704
    1.1049180327868853, 1255113.4937515096
    1.1852459016393442, 1034416.372625131
    1.2737704918032788, 886319.8712540427
    1.359016393442623, 779309.6059234422
    1.4540983606557378, 609961.3425615291
    1.4918032786885247, 502533.5248888627
    1.5098360655737708, 398203.214420678
    1.5196721311475412, 299601.3683763488
    1.5213114754098362, 195504.2893801002
    1.5262295081967214, 119588.7737216253
    1.5360655737704918, 65114.84416472715
    1.5475409836065577, 38316.26954492593
    1.5508196721311474, 25329.082065449114
    1.559016393442623, 18331.533361062826
    1.5688524590163935, 14713.897024069562
    1.5934426229508198, 11361.99263934681
    1.6196721311475408, 8773.786275376307
    1.6508196721311474, 6433.761487810411
    1.6754098360655738, 5032.800618487584
    1.7049180327868854, 3787.196421050415
    1.7327868852459019, 2886.942647505133
    1.7557377049180327, 2287.679320744783
    1.777049180327869, 1933.9100405245356];


%planks constant
h = 4.135e-15; %eV s
%speed of light
c = 3e8; %m/s

%absorption coef
alpha = data(:,2); %1/m

%wavelength
lambda = data(:,1); %microns

%energy
E = h*c./(lambda*1e-6);


figure(1)
plot(lambda,alpha)
xlabel('Wavelength (\mum)')
ylabel('\alpha (m^{-1})')
set(gca,'Yscale','log')

figure(2)
plot(E,alpha)
xlabel('Energy (eV)')
ylabel('\alpha (m^{-1})')
set(gca,'Yscale','log')


hb = 1.986e-25;

figure(3)
plot(lambda,alpha)
hAx(1)=axes;                                   % first axes
hAx(2)=axes('Position',hAx(1).Position, ...    % 2nd on top of first
    'color','none', ...                % so will show through
    'XAxisLocation','top', ...         % move XAxis to top
    'XDir','reverse');                 % and flip LR
hAx(1).XLim=[0.2 1.8];                         % first axes limits
eV=hb./hAx(1).XLim*6.242E18*1E6;               % limits of 2nd in eV from first
hAx(2).XLim=sort(eV);                          % match the limits
set(hAx,{'YScale'},{'log'},{'YLim'},{[1E3 1E8]});  % log y axes
fneV=@(lam) hb./lam*6.242E18*1E6;        % function for wavelength-->eV
fnWL=@(eV) hb./eV*6.242E18*1E6;          % and vice versa...convenience
hAx(2).XLim=hAx(1).XLim;                 % use same limits to match up to
inpt = [0.7:0.1:0.9 1:5];
outpt = fnWL([0.7:0.1:0.9 1:5]);
lamAx=fnWL([0.7:0.1:0.9 1:5]);           % wavelengths for the wanted eV values
% draw ticks at matching wavelength, label with eV value
set(hAx(2),'XTick',sort(lamAx),'XTickLabel',fneV(lamAx))
hAx(2).FontSize=8;                       % need a little more room
hAx(1).YTick=[];                         % remove superfluous axis ticks
%box on                                   % add the box around the figure

[dpb]

Because the axes aren't linearly related, you have to scale the nonlinear one manually with the proper x data values across the xlim range.

That's what the two anonymous functions are for -- convert one to the other. I forget which was the dominant linear here and don't have time to dig at the moment. I'll try to come back later after other commitments...

[Christopher Saltonstall]

If you comment out the "plot" command, the code works. However, when data is plotted, another set of axes seems to be overlayed.

[dpb]

Yes, that would happen.

The code was written to draw the axes then plot() into those two existing axes with the axes handle for each in turn.

The code sequence

figure
plot()

will create the axes object in which the data are plotted, then the next two lines will create the two axes which are customized as desired.

Don't create another axes; remove that plot() call and replace with

hL(1)=plot(hAx(1),...);
hL(2)=plot(hAx(2),...);

where the x,y data for each are the proper for each of the two axes,respectively, and in the proper x units (after the two axes are created, of course).

I'm not positive otomh right now whether you also need hold on for each when using the axes handle or not; it can't hurt and may be needed to keep autoscaling from changing limits on you...

[Christopher Saltonstall]

Thanks. I eventually figured that out. Now I'm trying to find out why the ticks dont quite line up. Specifically, 0.9 eV = ~1.4 um not 1.6 um as shown above.

[dpb]

I'll have to come back tonight but I'm thinking there may yet need to be another scaling that I left out of the above -- the above places the two endpoints correctly, but still linear between those based on the computed eV(lam) value. That returns the right numeric value, but that isn't the place it goes precisely on the nonlinear axes to actually correspond to the same location on the linear lambda axes.

It makes my head hurt to think about... :)

But, there's where the rub lies, scratch your head over the algebra and I'll look in again later.

[Christopher Saltonstall]

Solution

data = [0.35737704918032787, 92649850.48039015
    0.3819672131147541, 72475211.53588514
    0.41147540983606556, 55967196.301264346
    0.4229508196721312, 51126728.207795605
    0.4557377049180328, 47937263.23846833
    0.5016393442622952, 43803474.996365294
    0.539344262295082, 40544606.59504045
    0.5688524590163934, 34277014.89316264
    0.5950819672131147, 25461360.619704828
    0.6229508196721312, 16618285.037019765
    0.6655737704918033, 10037712.197086804
    0.7213114754098361, 6724703.556947659
    0.7950819672131147, 4169337.0144349397
    0.8836065573770491, 2722604.5977291227
    0.9557377049180328, 1997130.2156342946
    1.0344262295081967, 1523018.3932672704
    1.1049180327868853, 1255113.4937515096
    1.1852459016393442, 1034416.372625131
    1.2737704918032788, 886319.8712540427
    1.359016393442623, 779309.6059234422
    1.4540983606557378, 609961.3425615291
    1.4918032786885247, 502533.5248888627
    1.5098360655737708, 398203.214420678
    1.5196721311475412, 299601.3683763488
    1.5213114754098362, 195504.2893801002
    1.5262295081967214, 119588.7737216253
    1.5360655737704918, 65114.84416472715
    1.5475409836065577, 38316.26954492593
    1.5508196721311474, 25329.082065449114
    1.559016393442623, 18331.533361062826
    1.5688524590163935, 14713.897024069562
    1.5934426229508198, 11361.99263934681
    1.6196721311475408, 8773.786275376307
    1.6508196721311474, 6433.761487810411
    1.6754098360655738, 5032.800618487584
    1.7049180327868854, 3787.196421050415
    1.7327868852459019, 2886.942647505133
    1.7557377049180327, 2287.679320744783
    1.777049180327869, 1933.9100405245356];

fig = figure(5);
% setup bottom axis
ax = axes();
hold(ax);
ax.YAxis.Scale = 'log';
xlabel(ax, 'Wavelength ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14);
ylabel(ax, '$\alpha$ (m$^{-1}$)', 'Interpreter', 'latex', 'FontSize', 14);
% setup top axis
ax_top = axes(); % axis to appear at top
hold(ax_top);
ax_top.Position = ax.Position;
ax_top.YAxis.Visible = 'off';
ax_top.XAxisLocation = 'top';
% ax_top.XDir = 'reverse';
ax_top.Color = 'none';
xlabel(ax_top, 'Energy($e$V)', 'Interpreter', 'latex', 'FontSize', 14);
h = 4.135e-15; %eV s
c = 3e8; %m/s
lambda = linspace(0.2,1.8,9); %m
E = h*c./lambda*10^6; % 10^6 because lambda is in microns
% configure limits of bottom axis
ax.XLim = [lambda(1) lambda(end)];
ax.XTick = lambda;
ax.XAxis.TickLength = [0.015, 0.00];
ax.YAxis.TickLength = [0.02, 0.00];
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = linspace(0.2,1.8,17);
% configure limits and labels of top axis
y_ticks = [0.7 0.8 0.9 1 2 3 4 5];
lambda_y_tick = h*c./y_ticks*10^6;
ax_top.XLim = [lambda(1) lambda(end)];
ax_top.XTick = fliplr(lambda_y_tick);
ax_top.XTickLabel = compose('%1.1f', fliplr(y_ticks));
ax_top.XAxis.TickLength = [0.02, 0.00];
ax_top.XAxis.MinorTick = 'off';
plot(ax, data(:,1), data(:,2), 'r-', 'LineWidth', 3);