Diode Circuit Analysis

From Wikimization

function [vl Cdmax Cdmin Cdmed] = diodeMap(vi, Rs, Rl, verbose, fundfreq, Fs, precision)  %from Convex Optimization & Euclidean Distance Geometry
   mp.Digits(precision);  %Advanpix MCT
   k   = 1.380649e-23;    %boltzmann constant  % https://www.nist.gov/si-redefinition/kelvin-boltzmann-constant
   T   = 298;             %temperature kelvin  %Neudeck p.5                
   q   = 1.602176634e-19; %electron charge     %NIST
   is  = 76.9e-12;        %circuitlab diodeInduced
   eta = 1.45;            %circuitlab diodeInduced
   tau = 4.32e-6;         %transport time. max exponent (observed) e-4. Increasing capacitance observes more nonlinearity in vL FFT.
   Rd  = 0.042;           %diode internal resistance
   nVt = eta*k*T/q;       %[V] thermal voltage for particular diode characterized by eta
   isR = is*(Rs*Rl/(Rs + Rl) + Rd);

   dvd = 1000*randn(size(vi));  %perturbation of derivative proves that iteration converges to same median capacitance
   for i=1:5  %3 iterations will introduce randomness into mean impedance calculation
      Lw = exp((isR + vi*Rl/(Rl + Rs))/nVt)*isR.*(1 + tau*dvd/nVt)/nVt;
      vd = isR - nVt*Lambert_W(Lw) + vi*Rl/(Rl + Rs);
      dvd = [vd(2)-vd(end); vd(3:end)-vd(1:end-2); vd(1)-vd(end-1)]/2;  %approximate the derivative first time around, otherwise dvd becomes near 0. Assumes periodicity.   
      if verbose
         Cd = tau*is*exp(vd/nVt)/nVt;
         Cdmed = median(Cd);
         disp([num2eng(Cdmed,false,false,false,16) ' farad']);
      end
   end
   if verbose
      Cdmax = max(Cd);
      Cdmin = min(Cd);

      figure(20); plot(0:numel(Cd)-1, Cd);
      title('Cd over time'); xlabel('sample number'); ylabel('Farad'); set(get(gca,'ylabel'),'rotation',0);
      drawnow
   else
      Cd = tau*is*exp(vd/nVt)/nVt;
   end
   iC = Cd.*dvd*Fs;            %diffusion capacitor current. Normalization by Fs makes current independent of sample rate
   id = is*(exp(vd/nVt) - 1);  %ideal Shockley equation.  This is already independent of Fs.
   vl = (id + iC)*Rd + vd;

   if verbose
      figure(21); plot(0:numel(iC)-1, iC);
      title('Cd current'); xlabel('sample number'); ylabel('Amp'); set(get(gca,'ylabel'),'rotation',0);
      figure(22); plot(0:numel(id)-1, id);
      title('diode current'); xlabel('sample number'); ylabel('Amp'); set(get(gca,'ylabel'),'rotation',0);
      idx = find(iC ~= 0);
      disp(['mean   capacitor  impedance  = ' num2str(orosumvec(vd(idx)./iC(idx),1)/numel(idx),'%3.4e')])
%       disp(['mean capacitor |impedance| = ' num2str(orosumvec(abs(vd(idx)./iC(idx)),1)/numel(idx),'%3.4e')])
      t = median(vd./iC);
      if t < 0, spc=[]; else spc=' '; end
      disp(['median capacitor  impedance  =' spc num2str(t)])  %median impedance is always far smaller than mean despite mp.Digits precision
      medcapimp = median(abs(vd./iC));  
      disp(['median capacitor |impedance| = ' num2str(medcapimp,'%3.4e')])
      disp(['empirical constant c = ' num2str(medcapimp*(2*pi*fundfreq*Cdmed),'%3.15e')])

      figure(23);
      histogram(Cd,50,'Normalization','probability');                    %for book illustration
   %    histogram(Cd,round(3.25*48000));  %for counting smallest diffusion capacitors
      title('histogram C_d'); xlabel('bin'); ylabel('relative count'); %set(get(gca,'ylabel'),'rotation',0);
      drawnow
   end
end

% https://www.mathworks.com/matlabcentral/fileexchange/43419-the-lambert-w-function
function w = Lambert_W(x,branch)
% Lambert_W  Functional inverse of x = w*exp(w).
% w = Lambert_W(x), same as Lambert_W(x,0)
% w = Lambert_W(x,0)  Primary or upper branch, W_0(x)
% w = Lambert_W(x,-1)  Lower branch, W_{-1}(x)
%
% See: http://blogs.mathworks.com/cleve/2013/09/02/the-lambert-w-function/

   % Effective starting guess
   if nargin < 2 || branch ~= -1
      w = ones(size(x));     % Start above -1
   else  
      w = -2*ones(size(x));  % Start below -1
   end
   v = inf*w;

   % Halley's method
   c = 0;
   limit = 100;
   while any(abs(w - v)./abs(w) > 1e-15) && c < limit  %changed magic tolerance. Was 1e-8. -JonD
      v = w;
      e = exp(w);
      f = w.*e - x;  % Iterate, to make this quantity zero
      w = w - f./((e.*(w+1) - (w+2).*f./(2*w + 2)));
      c = c + 1;
   end
   if c >= limit
      disp('%%%%%%%%%%%%%%%%%%%%%% Warning: Lambert limit reached %%%%%%%%%%%%%%%%%%%%%')
   end
end