%Prof. E.R. Brown, 16 January 2008, composed for ECE215, "Solid State
%Engineering for Electronics."
%This solves for the locus of points at which an electron in an ionized
%hydrogen molecule provides an attractive effect on the protons, with and
%without proton-proton repulsion included.  Only Coulomb effects are
%considered, whicn lead to an implicit hyperbolic equation that is solved by
%by iteration using MATLAB in an x-y coordinate system with (x,y)=(0,0)
%at the midpoint between the two protons and avoiding the y=0,x = +/-a/2 
%points to avoid the Coulombic singularity.  It solves under two conditions: 
%(1) ignoring the proton-proton repulsion (output plot in blue), 
%and (2) including the proton-proton repulsion (output plot in red).
%It only solves in the upper-right quadrant, 
%using symmetry to construct the entire curve in the other three quadrants.
%Input variables: a is the proton-proton separation;
%xstart should be 0.0, xstop should be a or greater;
%make sure ystart is > 0 to avoid the Coulomb singularity at the
%proton locations; ystop should be somewhat greater than a/2.
%a typical input command from the workspace is:
%hydrogen(0.0,2.0,0.01,2.0,2.0)


function f = hydrogen(xstart,xstop,ystart,ystop,a);


dx=(xstop-xstart)./1000
dy=(ystop-ystart)./1000
%first do the analysis for electron-proton interaction only
val1=0
for iy=1:1000
    y=ystart+dy*(iy-1);
    for ix=1:1000
        x=xstart+dx*(ix-1);
val2 = (x+a/2)./[(x+a/2).^2 + y.^2].^1.5 - (x-a/2)./[(x-a/2).^2 + y.^2].^1.5;
 if (val1>0)&(val2<0)
           locus1(iy,:)=[x,y];
           locus2(iy,:)=[-x,y];
           locus3(iy,:)=[-x,-y];
           locus4(iy,:)=[x,-y];
 end
        val1=val2;
    end

end
%now do the computation fo electron-proton and proton-proton interaction
val1=0
for iy=1:1000
    y=ystart+dy*(iy-1);
    for ix=1:1000
        x=xstart+dx*(ix-1);
val2 = (x+a/2)./[(x+a/2).^2 + y.^2].^1.5 - (x-a/2)./[(x-a/2).^2 + y.^2].^1.5 - 1./a.^2;
 if (val1>0)&(val2<0)
           locus5(iy,:)=[x,y];
           locus6(iy,:)=[-x,y];
           locus7(iy,:)=[-x,-y];
           locus8(iy,:)=[x,-y];
 end
        val1=val2;
    end

end
    figure(1)
    plot(locus1(:,1),locus1(:,2),'b');
    hold on
    plot(locus2(:,1),locus2(:,2),'b');
    plot(locus3(:,1),locus3(:,2),'b');
    plot(locus4(:,1),locus4(:,2),'b');
    plot(locus5(:,1),locus5(:,2),'r');
    plot(locus6(:,1),locus6(:,2),'r');
    plot(locus7(:,1),locus7(:,2),'r');
    plot(locus8(:,1),locus8(:,2),'r');
    plot(a/2,0,'+k','MarkerSize',12);
    plot(a/2,0,'ok','MarkerSize',12);
    plot(-a/2,0,'+k','MarkerSize',12);
    plot(-a/2,0,'ok','MarkerSize',12);
    
    clear all
%fprintf(sanalyzer,'FA:?') %ask for start frequency
%start_freq=str2num(fscanf(sanalyzer))
%subplot(1,1,1)
 %   figure(nloop+5),plot(time,SNRpower,'*')
 %   xlabel('Time [ns]')
 %   ylabel('Channel SNR ')
    
  %  figure(nloop+6),plot(time,SNRcross,'*')
  %  xlabel('Time [ns]')
  %  ylabel('Cross Channel SNR ')