如何用颜色强度在3D空间表示p[r0,gamma,a]大小？R实现咨询

完全可行，rgl包的plot3d或points3d函数支持通过col参数指定每个点的颜色，只需将三维坐标和对应的p值转换为向量形式，再将p值映射到颜色空间即可。以下是具体实现步骤和代码：

步骤说明 #

1. 处理p的数值：你的p数组存储的是复数形式的声压，需先计算其幅值（绝对值）作为颜色对应的量化值，因为颜色只能映射实数。
2. 展开三维坐标：将r0、gamma、a的所有组合转换为三个一维向量，与展开后的p幅值向量一一对应。
3. 颜色映射：将p的幅值归一化后，映射到渐变颜色（比如从蓝色到红色的热图色系）。
4. 绘制三维散点图：使用plot3d，传入坐标向量和颜色向量。

修正并补充后的完整代码 #

library(rgl)
library(scales) # 用于归一化数值

# 修正p数组维度：循环中a只到10，所以dim设为c(20,36,10)
p = array(0, dim=c(20,36,10))
q = c(1:20)

k = 4187         # phase constant [rad/m] = (2*pi)/lambda =(2*3.14)/(1500/1000000)
f = 1000000   # frequency [Hz]
rou = 1.06      # density

for (r0 in 1:20)     # Depth [*5 mm]
{
  for (gamma in 1:36)
  {
    for (a_idx in 1:10)  # 用a_idx避免覆盖循环变量a
    {
      a = a_idx*0.0005
      Jn_z = 0
      for (m in 1:100)
      {
        z = k*a*sin(gamma*10 * pi/180) # 注意：sin函数默认弧度，这里gamma*10是角度，需转弧度
        P_Jn = (((-1)^m)/(factorial(m)*factorial(m+1)))*(z/2)^(2*m+1)
        Jn_z = Jn_z + P_Jn
      }
      Imp<-complex(real=0,imaginary=2*pi*f*rou)
      Phase<-complex(real=0,imaginary=2*pi*f - k*(r0*5e-3)) # r0是*5mm，转米
      
      # Sound Pressure 
      p[r0,gamma,a_idx] = Imp*(a^2*5*exp(Phase))/( (r0*5e-3)*(Jn_z/(k*a*sin(gamma*10 * pi/180))) )
    }
  }
}

# ---------------------- 绘图部分 ----------------------
# 1. 计算p的幅值
p_amp = Mod(p)

# 2. 生成所有坐标组合
# 转换为实际物理值：r0是深度（米），gamma是角度（度），a是膜片半径（米）
coords = expand.grid(
  r0 = (1:20)*5e-3,
  gamma = (1:36)*10,
  a = (1:10)*0.0005
)

# 3. 展开p_amp为向量，与coords顺序对应
p_amp_vec = as.vector(p_amp)

# 4. 颜色映射：归一化幅值到0-1，生成热图色系
color_palette = colorRampPalette(c("blue", "green", "yellow", "red"))(100)
normalized_amp = rescale(p_amp_vec, to = c(1, 100))
point_colors = color_palette[round(normalized_amp)]

# 5. 绘制三维散点图
plot3d(
  x = coords$r0, y = coords$gamma, z = coords$a,
  col = point_colors,
  size = 2, # 调整点大小
  xlab = "Depth (m)", ylab = "Gamma (deg)", zlab = "Diaphragm Radius (m)",
  main = "Sound Pressure Amplitude in 3D Space"
)

# 添加颜色图例
legend3d("topright", legend = c("Low", "Medium", "High"),
         col = color_palette[c(1, 50, 100)], pch = 16, cex = 0.8)

关键细节说明 #

• 修正了原代码中的几个潜在问题：
  • p数组维度与循环范围匹配（从100改为10）
  • sin函数使用弧度计算（原代码直接用角度值传入，结果错误）
  • r0转换为实际米单位参与计算
  • 用a_idx作为循环变量，避免覆盖a的物理值
• 颜色映射使用colorRampPalette生成渐变色系，你可以根据需求替换为其他色系（如rainbow、heat.colors）
• scales::rescale用于将幅值归一化到颜色索引的范围，确保颜色渐变均匀
• legend3d用于添加颜色图例，方便解读颜色对应的幅值大小

内容的提问来源于stack exchange，提问作者okadasan