docs: def densu(alt, dlb, tinf, tlb, xm, alpha, tz, zlb, s2, mn1, zn1, tn1, tgn1):
      '''
  /*      Calculate Temperature and Density Profiles for MSIS models
   *      New lower thermo polynomial
   */
   tz, zn1, tn1, and tgn1 are simulated pointers
   '''
      rgas = 831.4
      #rgas = 831.44621    #maybe make this a global constant?
      densu_temp = 1.0
  
      xs = [0.0]*5
      ys = [0.0]*5
      y2out = [0.0]*5
  
      #/* joining altitudes of Bates and spline */
      za=zn1[0];
      if (alt>za):
          z=alt;
      else:
          z=za;
  
      #/* geopotential altitude difference from ZLB */
      zg2 = zeta(z, zlb);
      
      #/* Bates temperature */
      tt = tinf - (tinf - tlb) * exp(-s2*zg2);
      ta = tt;
      tz[0] = tt
      densu_temp = tz[0]
  
      if (alt<za):
          #/* calculate temperature below ZA
          # * temperature gradient at ZA from Bates profile */
          dta = (tinf - ta) * s2 * pow(((re[0]+zlb)/(re[0]+za)),2.0);
          tgn1[0]=dta;
          tn1[0]=ta;
          if (alt>zn1[mn1-1]):
              z=alt;
          else:
              z=zn1[mn1-1];
          mn=mn1;
          z1=zn1[0];
          z2=zn1[mn-1];
          t1=tn1[0];
          t2=tn1[mn-1];
          #/* geopotental difference from z1 */
          zg = zeta (z, z1);
          zgdif = zeta(z2, z1);
          #/* set up spline nodes */
          for k in range(mn):
              xs[k] = zeta(zn1[k], z1) / zgdif;
              ys[k] = 1.0 / tn1[k];
          
          #/* end node derivatives */
          yd1 = -tgn1[0] / (t1*t1) * zgdif;
          yd2 = -tgn1[1] / (t2*t2) * zgdif * pow(((re[0]+z2)/(re[0]+z1)),2.0);
          #/* calculate spline coefficients */
          spline (xs, ys, mn, yd1, yd2, y2out);
          x = zg / zgdif;
          y = [0.0]
          splint (xs, ys, y2out, mn, x, y);
          #/* temperature at altitude */
          tz[0] = 1.0 / y[0];
          densu_temp = tz[0];
  
      if (xm==0):
          return densu_temp;
  
      #/* calculate density above za */
      glb = gsurf[0] / pow((1.0 + zlb/re[0]),2.0);
      gamma = xm * glb / (s2 * rgas * tinf);
      expl = exp(-s2 * gamma * zg2);
      if (expl>50.0): # pragma: no cover
          expl=50.0;
      if (tt<=0): # pragma: no cover
          expl=50.0;
  
      #/* density at altitude */
      densa = dlb * pow((tlb/tt),((1.0+alpha+gamma))) * expl;
      densu_temp=densa;
      if (alt>=za):
          return densu_temp;
  
      #/* calculate density below za */
      glb = gsurf[0] / pow((1.0 + z1/re[0]),2.0);
      gamm = xm * glb * zgdif / rgas;
  
      #/* integrate spline temperatures */
      yi = [0]
      splini (xs, ys, y2out, mn, x, yi);
      expl = gamm * yi[0];
      if (expl>50.0): # pragma: no cover
          expl=50.0;
      if (tz[0]<=0): # pragma: no cover
          expl=50.0;
  
      #/* density at altitude */
      densu_temp = densu_temp * pow ((t1 / tz[0]),(1.0 + alpha)) * exp(-expl);
      return densu_temp;