Code doesn't do calculation

[larjy]

Overview: I'm creating a code to simulate quantum tunneling current.

I'm getting the code to output data on a spreadsheet file, and it all works except the following part, which are just giving values of 0 all the way down for the entire column.

j_tl=((eV*hbar*k)/(m/c*c))*(1 - R_lr - T_rl);

When I remove the ((eV*hbar*k)/(m/c*c)) part of it, it gives out values for the (1 - R_lr - T_rl) so I presume the problem is in the former. But I'm struggling to think what it could be. The same is happening for j_tr.

I would really appreciate some help. I feel it may be something trivial that I'm just missing.

Full code here

#include <iostream>
#include <math.h>
#include <fstream>
#include <stdio.h>

using namespace std;

int main()

{
    double k,q,p,y;
    double E; //Electron energy (V)
    double V=0.3; //Bias (V)
    double U=1.5; //Barrier Height (V)
    double d=3.0e-10; //Barrier width (m)

    const double m=0.510998928e6; //Electron mass (eV/c^2)
    const double hbar=6.58211928e-16; //Reduced Planck constant (eV.s)
    const double c=299792458.0; //Speed of light (m/s)
	const double eV=1.60217657e-19; //Electron charge

    double R_lr, R_rl; //Reflection coefficient
    double T_lr, T_rl; //Transmission coefficient

	double j_tl; //Total Current left region
	double j_tr; //Total Current right region

    FILE *STMdata;
    STMdata=fopen("STMdata_table.xls","w");
    fprintf(STMdata,"Energy \t Reflection (L-R) \t Transmission (L-R) \t Reflection (R-L) \t Transmission (R-L) \t Total Current Left \t Total Current Right \n");

    for(float E=0.5;E<10.0;E+=0.01)
    {

    if (E>U+V)
    {
        k=((sqrt(2*(m/(c*c))*(E-V)))/hbar);
        q=((sqrt(2*(m/(c*c))*(E-U-V)))/hbar);
        p=((sqrt(2*(m/(c*c))*E))/hbar);

        R_lr=(q*q*(k-p)*(k-p)*cos(q*d)*cos(q*d)+(q*q-p*k)*(q*q-p*k)*sin(q*d)*sin(q*d))/
             (q*q*(k+p)*(k+p)*cos(q*d)*cos(q*d)+(q*q+p*k)*(q*q+p*k)*sin(q*d)*sin(q*d));

        T_lr=(4*k*k*q*q)/
             (q*q*(k+p)*(k+p)*cos(q*d)*cos(q*d)+(q*q+p*k)*(q*q+p*k)*sin(q*d)*sin(q*d));

		R_rl=R_lr;

		T_rl=((p*p)/(k*k))*T_lr;

		//Total Current

		j_tl=((eV*hbar*k)/(m/c*c))*(1 - R_lr - T_rl);

		j_tr=((eV*hbar*p)/(m/c*c))*(T_lr + R_rl - 1);
    }

    else if (V<E<U+V)
    {
	    k=((sqrt(2*(m/(c*c))*(E-V)))/hbar);
        y=((sqrt(2*(m/(c*c))*(V+U-E)))/hbar);
        p=((sqrt(2*(m/(c*c))*E))/hbar);

        R_lr=(y*y*(k-p)*(k-p)*cosh(y*d)*cosh(y*d)+(-(y*y)-p*k)*(-(y*y)-p*k)*sinh(y*d)*sinh(y*d))/
             (y*y*(k+p)*(k+p)*cosh(y*d)*cosh(y*d)+(-(y*y)+p*k)*(-(y*y)+p*k)*sinh(y*d)*sinh(y*d));

        T_lr=(4*k*k*y*y)/
             (y*y*(k+p)*(k+p)*cosh(y*d)*cosh(y*d)+(-(y*y)+p*k)*(-(y*y)+p*k)*sinh(y*d)*sinh(y*d));

		R_rl=R_lr;

		T_rl=((p*p)/(k*k))*T_lr;

		//Total Current

		j_tl=((eV*hbar*k)/(m/c*c))*(1 - R_lr - T_rl);

		j_tr=((eV*hbar*p)/(m/c*c))*(T_lr + R_rl - 1);

    }


    fprintf(STMdata,"%f \t %f \t %f \t %f \t %f \t %f \t %f \n",E,R_lr,T_lr,R_rl, T_rl,j_tl,j_tr);

    }

    fclose(STMdata);


    return 0;
}

[MiiNiPaa]

if (V<E<U+V)
How is that works:

V < E < (U+V)
//        ↓(result we will call x)
(V < E) < x
// ↓(result will be bool which converts to 0 or 1)
(0 or 1) < x
//↑↑resulting expression 

You cannot chain comparsion operators. You should write if (V < E && E < U + V)

[larjy]

Thank you for your reply.

I didn't realise that V<E<U+V would pose a problem, it had worked fine for giving all the other values.

I changed it to your suggestion, problem still persists, it's made no difference.

[MiiNiPaa]

Another probably undesired behavior: m/c*c is calculated as (m/c)*c and equals to m

I will not check all of your code for similar errors, do it yourself.

Another question: lines 41 and 45: you are calculating R_lr and T_lr using a normal cos.
lines 64 and 67: you are calculating same variables using hyperbolic cos. Is it supposed to be that way?

[larjy]

Thanks for your time.

Made that change too, it hasn't done anything.

The fundamental code works, 4 of the 6 data columns have the right output. Even with the trivial language errors, I should get some (all be it incorrect) values for the remaining 2 columns using that j_tl and j_tr calculation but I instead am getting a list of 0's.

If I apply the calculation manually on excel using the data from the 4 columns to get columns 5 and 6 I get the right values. The calculation is fairly straight forward but for whatever reason it's just not coming together.

Yes it's supposed to be hyperbolic cos.

Anyway I really appreciate your time, thanks for the suggestions.

[MiiNiPaa]

Found it: %f format specifier prints value in decimal format with 6 digit precision. And your values are very small, so first 6 digits are 0. Use %e instead.

You can also look into %g to avoid need to chack possible output range: http://en.cppreference.com/w/cpp/io/c/fprintf

[larjy]

I can't thank you enough!

That's fixed it. Brilliant, I really appreciate your help.