
macro find.fwhm.pgauss
  if ($0 != 1)
    echo "USAGE: find.fwhm"
    break
  end
  
  $dz = 0.01
  create z 0 50 $dz

  # Gaussian taylor expansion
  # f = 1 / (1 + z + z^2/2 + z^3/6)
  # 1 + z + z^2/2 + z^3/6 = 2
  # f = z + z^2/2 + z^3/6 - 2, find where f == 0.0
  
  set f = z + 0.5*z^2 + (1/6.0)*z^3 - 2.0
  set dfdz = 1 + z + 0.5*z^2

  lim -n 0 z f; clear; box; plot z f
  lim -n 1 z dfdz; clear; box; plot z dfdz  

  $nZ0 = 0
  $nZ1 = 5 / $dz

  $Zg = 0.5*(z[$nZ0] + z[$nZ1])  
  $nZg = int($Zg / $dz)
  $dZ = 1.0

  for i 0 10
    $Fg = f[$nZg]
    $dFdz_g = dfdz[$nZg]

    $dZ = $Fg / $dFdz_g

    echo $Zg $Fg $dZ $dFdz_g

    $Zg -= $dZ
    $Zg = max ($Zg , 0)
    $nZg = int($Zg / $dz)
  end 
  echo $Zg $Fg $dZ $dFdz_g
  $Zhm = $Zg
  $FWHM = 2*sqrt(2*$Zg)

  echo $Zhm : $FWHM
end

macro find.fwhm.rgauss
  if ($0 != 2)
    echo "USAGE: find.fwhm (K)"
    break
  end
  
  $K = $1

  if ($K == 0.0)
    $Zhm = (sqrt(5) - 1.0) / 2.0
    $FWHM = 2*sqrt(2*$Zhm)
    echo $K : $Zhm : $FWHM
    return
  end

  $dz = 0.01
  create z 0 50 $dz

  # set f = 1.0 / (1.0 + $K*z + z^1.667)
  # f = $K*z + z^1.667 - 1.0, find where f == 0.0
  
  set f = z + z^$K - 1.0
  set dfdz = ln(z) * z^$K + 1.0

  lim -n 0 z f; clear; box; plot z f
  lim -n 1 z dfdz; clear; box; plot z dfdz  

  $nZ0 = 0
  $nZ1 = 5 / $dz

  $Zg = 0.5*(z[$nZ0] + z[$nZ1])  
  $nZg = int($Zg / $dz)
  $dZ = 1.0

  for i 0 10
    $Fg = f[$nZg]
    $dFdz_g = dfdz[$nZg]

    $dZ = $Fg / $dFdz_g

    echo $nZg $Zg $Fg $dZ $dFdz_g

    $Zg -= $dZ
    $Zg = max ($Zg , 0)
    $nZg = int($Zg / $dz)
  end 
  # echo $Zg $Fg $dZ $dFdz_g
  $Zhm = $Zg
  $FWHM = 2*sqrt(2*$Zg)

  echo $K : $Zhm : $FWHM
end

macro find.fwhm.ps1v1
  if ($0 != 2)
    echo "USAGE: find.fwhm (K)"
    break
  end
  
  $K = $1

  if ($K == 0.0)
    $Zhm = 1.0
    $FWHM = 2*sqrt(2)
    echo $K : $Zhm : $FWHM
    return
  end

  $dz = 0.01
  create z 0 50 $dz

  # set f = 1.0 / (1.0 + $K*z + z^1.667)
  # f = $K*z + z^1.667 - 1.0, find where f == 0.0
  
  set f = $K*z + z^1.667 - 1.0
  set dfdz = $K + z^0.667

  #lim -n 0 z f; clear; box; plot z f
  #lim -n 1 z dfdz; clear; box; plot z dfdz  

  $nZ0 = 0
  $nZ1 = 5 / $dz

  $Zg = 0.5*(z[$nZ0] + z[$nZ1])  
  $nZg = int($Zg / $dz)
  $dZ = 1.0

  for i 0 10
    $Fg = f[$nZg]
    $dFdz_g = dfdz[$nZg]

    $dZ = $Fg / $dFdz_g

    # echo $Zg $Fg $dZ $dFdz_g

    $Zg -= $dZ
    $Zg = max ($Zg , 0)
    $nZg = int($Zg / $dz)
  end 
  # echo $Zg $Fg $dZ $dFdz_g
  $Zhm = $Zg
  $FWHM = 2*sqrt(2*$Zg)

  echo $K : $Zhm : $FWHM
end

# qgauss is like ps1_v1 with z^2.25
macro find.fwhm.qgauss
  if ($0 != 2)
    echo "USAGE: find.qgauss (K)"
    break
  end
  
  $K = $1

  if ($K == 0.0)
    $Zhm = 1.0
    $FWHM = 2*sqrt(2)
    echo $K : $Zhm : $FWHM
    return
  end

  $dz = 0.01
  create z 0 50 $dz

  # set f = 1.0 / (1.0 + $K*z + z^2.25)
  # f = $K*z + z^2.25 - 1.0, find where f == 0.0
  
  set f = $K*z + z^2.25 - 1.0
  set dfdz = $K + z^1.25

  #lim -n 0 z f; clear; box; plot z f
  #lim -n 1 z dfdz; clear; box; plot z dfdz  

  $nZ0 = 0
  $nZ1 = 5 / $dz

  $Zg = 0.5*(z[$nZ0] + z[$nZ1])  
  $nZg = int($Zg / $dz)
  $dZ = 1.0

  for i 0 10
    $Fg = f[$nZg]
    $dFdz_g = dfdz[$nZg]

    $dZ = $Fg / $dFdz_g

    # echo $Zg $Fg $dZ $dFdz_g

    $Zg -= $dZ
    $Zg = max ($Zg , 0)
    $nZg = int($Zg / $dz)
  end 
  # echo $Zg $Fg $dZ $dFdz_g
  $Zhm = $Zg
  $FWHM = 2*sqrt(2*$Zg)

  echo $K : $Zhm : $FWHM
end

macro fwhm.trend

  delete fwhm_v k_v

  for k 0 20 1.0
    find.fwhm.qgauss $k
    concat $k k_v
    concat $FWHM fwhm_v
  end
end

