docs: def densm(alt, d0, xm, tz, mn3, zn3, tn3, tgn3, mn2, zn2, tn2, tgn2):
      '''
  /*      Calculate Temperature and Density Profiles for lower atmos.  */
  '''
      xs = [0.0]*10
      ys = [0.0]*10
      y2out = [0.0]*10
      rgas = 831.4
      #rgas = 831.44621    #maybe make this a global constant?
      densm_tmp=d0
      if (alt>zn2[0]): # pragma: no cover
          if(xm==0.0):
              return tz[0]
          else:
              return d0
  
      #/* STRATOSPHERE/MESOSPHERE TEMPERATURE */
      if (alt>zn2[mn2-1]):
          z=alt
      else:
          z=zn2[mn2-1];
      mn=mn2;
      z1=zn2[0];
      z2=zn2[mn-1];
      t1=tn2[0];
      t2=tn2[mn-1];
      zg = zeta(z, z1);
      zgdif = zeta(z2, z1);
  
      #/* set up spline nodes */
      for k in range(mn):
          xs[k]=zeta(zn2[k],z1)/zgdif;
          ys[k]=1.0 / tn2[k];
      yd1=-tgn2[0] / (t1*t1) * zgdif;
      yd2=-tgn2[1] / (t2*t2) * zgdif * (pow(((re[0]+z2)/(re[0]+z1)),2.0));
  
      #/* calculate spline coefficients */
      spline (xs, ys, mn, yd1, yd2, y2out);   #No need to change this
      x = zg/zgdif;
      y = [0.0]
      splint (xs, ys, y2out, mn, x, y);
  
      #/* temperature at altitude */
      tz[0] = 1.0 / y[0];
      if (xm!=0.0):
          #/* calaculate stratosphere / mesospehere density */
          glb = gsurf[0] / (pow((1.0 + z1/re[0]),2.0));
          gamm = xm * glb * zgdif / rgas;
  
          #/* Integrate temperature profile */
          yi = [0.0]
          splini(xs, ys, y2out, mn, x, yi);
          expl=gamm*yi[0];
          if (expl>50.0): # pragma: no cover
              expl=50.0
  
          #/* Density at altitude */
          densm_tmp = densm_tmp * (t1 / tz[0]) * exp(-expl);
  
      if (alt>zn3[0]):
          if (xm==0.0):
              return tz[0]
          else:
              return densm_tmp
  
      #/* troposhere / stratosphere temperature */
      z = alt;
      mn = mn3;
      z1=zn3[0];
      z2=zn3[mn-1];
      t1=tn3[0];
      t2=tn3[mn-1];
      zg=zeta(z,z1);
      zgdif=zeta(z2,z1);
  
  
  
      #/* set up spline nodes */
      for k in range(mn):
          xs[k] = zeta(zn3[k],z1) / zgdif;
          ys[k] = 1.0 / tn3[k];
      
      yd1=-tgn3[0] / (t1*t1) * zgdif;
      yd2=-tgn3[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];
      if (xm!=0.0):
          #/* calaculate tropospheric / stratosphere density */
          glb = gsurf[0] / (pow((1.0 + z1/re[0]),2.0));
          gamm = xm * glb * zgdif / rgas;
  
          #/* Integrate temperature profile */
          yi = [0.0]
          splini(xs, ys, y2out, mn, x, yi);
          expl=gamm*yi[0];
          if (expl>50.0): # pragma: no cover
              expl=50.0;
  
          #/* Density at altitude */
          densm_tmp = densm_tmp * (t1 / tz[0]) * exp(-expl);
      
      if (xm==0.0):
          return tz[0];
      else:
          return densm_tmp;