API tools faq
paste
Guest User

Untitled

a guest
Sep 30th, 2025
19
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 7.43 KB | Science | 0 0
raw download clone embed print report
  1. -
  2. - using LinearAlgebra, QuadGK, Distributions, StaticArrays, StatsBase, Parameters, Plots, Cubature
  3. -
  4. - function dich(circuit)
  5. 0 f = 7.5e6
  6. - #circuit = (C1=100e-12, C2=100e-12, L1=40e-9, L2=50e-9, Lc=1e-9, V1=100, Vg=20)
  7. - T = 1e-3
  8. - s = 2π / T + 2π * f * 1im
  9. -
  10. 0 @unpack C1, C2, L1, L2, Lc, V1, Vg = circuit
  11. -
  12. - Γg₁, ΓVth₁, θ₁ = Gamma(182.3, 0.0303), Gamma(172.27, 0.04218), -7.5
  13. - Γg₂, ΓVth₂, θ₂ = Gamma(182.3, 0.0303), Gamma(172.27, 0.04218), -7.5
  14. -
  15. - # Define vectors and matrices
  16. -
  17. 0 A = @SMatrix [1/(s*L1)+s*C1 0 -s*C1 -1/(s*L1) 0; 0 1/(s*L2)+s*C2 -s*C2 -1/(s*L2) 0; -s*C1 -s*C2 s*(C1+C2) 0 1; -1/(s*L1) -1/(s*L2) 0 1/(s*Lc)+1/(s*L1)+1/(s*L2) 0; 0 0 1 0 0]
  18. 0 B(g) = @SMatrix [-g[1] 0 0 0 0; 0 -g[2] 0 0 0; g[1] g[2] 0 0 0; 0 0 0 0 0; 0 0 0 0 0]
  19. - MI = @SMatrix [1 0 0 0 0; 0 1 0 0 0]
  20. 0 mD(g) = @SMatrix [g[1]*Vg/s; g[2]*Vg/s; -(g[1] + g[2])*Vg/s; 0; 0] # try adding views here and . dots
  21. 0 MD(g) = @SMatrix [-g[1]/s 0; 0 -g[2]/s; g[1]/s g[2]/s; 0 0; 0 0]
  22. - #Iinit = @SMatrix [-C1 * (V1 - Vg + mean(ΓVth₁)); -C2 * (V1 - Vg + mean(ΓVth₂)); (C1 + C2) * (V1 - Vg) + C1 * mean(ΓVth₁) + C2 * mean(ΓVth₂); 0; V1 / s]
  23. 0 Iinit = @SMatrix [-C1 * (V1 - Vg + mean(Gamma(172.27, 0.04218))); -C2 * (V1 - Vg + mean(Gamma(172.27, 0.04218))); (C1 + C2) * (V1 - Vg) + C1 * mean(Gamma(172.27, 0.04218)) + C2 * mean(Gamma(172.27, 0.04218)); 0; V1 / s]
  24. -
  25. - # Define distributions
  26. -
  27. 0 frankcopula_density(u₁, u₂, θ) = θ * (1 - exp(-θ)) * exp(-θ * (u₁ + u₂)) * ((1 - exp(-θ)) - (1 - exp(-θ * u₁)) * (1 - exp(-θ * u₂)))^-2
  28. -
  29. - # Finv_Vth1(u) = quantile(ΓVth₁, u)
  30. - # Finv_g1(u) = quantile(Γg₁, u)
  31. - # Finv_Vth2(u) = quantile(ΓVth₂, u)
  32. - # Finv_g2(u) = quantile(Γg₂, u)
  33. - # this halves allocations, why?
  34. 0 @inline Finv_Vth1(u) = quantile(Gamma(172.27, 0.04218), u)
  35. 0 @inline Finv_g1(u) = quantile(Gamma(182.3, 0.0303), u)
  36. 0 @inline Finv_Vth2(u) = quantile(Gamma(172.27, 0.04218), u)
  37. 0 @inline Finv_g2(u) = quantile(Gamma(182.3, 0.0303), u)
  38. -
  39. 0 f_Vth₁_g₁(v, g) = frankcopula_density(cdf(ΓVth₁, v), cdf(Γg₁, g), θ₁) * pdf(ΓVth₁, v) * pdf(Γg₁, g)
  40. 0 f_Vth₂_g₂(v, g) = frankcopula_density(cdf(ΓVth₂, v), cdf(Γg₂, g), θ₂) * pdf(ΓVth₂, v) * pdf(Γg₂, g)
  41. - f_Vth_g(v, g) = f_Vth₁_g₁(v[1], g[1]) * f_Vth₂_g₂(v[2], g[2])
  42. -
  43. - # Calculate matrices
  44. -
  45. 0 K1(g) = MI * (A * inv(A - B(g)) * (mD(g) + Iinit) - Iinit)
  46. 0 K2(g) = MI * A * (inv(A - B(g))) * MD(g)
  47. -
  48. - # Define h-tilde and h-tilde^-1
  49. -
  50. 0 h(Vth, g) = K1(g) + K2(g) * Vth
  51. -
  52. 0 function hinvg(Ich)
  53. 0 γ = real(1/s) / imag(1/s)
  54. 0 Ĩch = @SMatrix [Ich[1]; Ich[2]; -Ich[1]-Ich[2]; 0; 0]
  55. 0 u0 = real(MI * Ĩch) - γ * imag(MI * Ĩch)
  56. 0 invA_I = inv(A) * (Ĩch + Iinit)
  57. 0 U_1 = real(MI * B(@SMatrix [1 0]) * invA_I) - γ * imag(MI * B(@SMatrix [1 0]) * invA_I)
  58. 0 U_2 = real(MI * B(@SMatrix [0 1]) * invA_I) - γ * imag(MI * B(@SMatrix [0 1]) * invA_I)
  59. - U = [U_1 U_2]
  60. 0 g = inv(U) * u0
  61. 0 return g
  62. - end
  63. -
  64. 0 hinvVth(Ich, g) = real(inv(K2(g)) * (Ich - K1(g))) # real() used to avoid numbers become complex by numerical calculations
  65. -
  66. 0 function hinv(Ich)
  67. - g = hinvg(Ich)
  68. - Vth = hinvVth(Ich, g)
  69. - return (Vth, g)
  70. - end
  71. -
  72. - # Define Jacobian
  73. -
  74. 0 function dhinv(Vth, g)
  75. 0 dMDdg1 = @SMatrix [-1/s 0; 0 0; 1/s 0; 0 0; 0 0]
  76. 0 dMDdg2 = @SMatrix [0 0; 0 -1/s; 0 1/s; 0 0; 0 0]
  77. 0 dmDdg1 = @SMatrix [Vg/s; 0; -Vg/s; 0; 0]
  78. 0 dmDdg2 = @SMatrix [0; Vg/s; -Vg/s; 0; 0]
  79. 0 invA_B = inv(A - B(g))
  80. 0 invA_B_M = invA_B * (mD(g) + MD(g) * Vth + Iinit)
  81. 0 dhdg1 = MI * A * (invA_B * (B(@SMatrix [1 0])) * invA_B_M + invA_B * (dmDdg1 + dMDdg1 * Vth))
  82. 0 dhdg2 = MI * A * (invA_B * (B(@SMatrix [0 1])) * invA_B_M + invA_B * (dmDdg2 + dMDdg2 * Vth))
  83. 0 Jh = vcat(hcat(real(K2(g)), real(dhdg1), real(dhdg2)), hcat(imag(K2(g)), imag(dhdg1), imag(dhdg2)))
  84. 0 return inv(Jh)
  85. - end
  86. -
  87. 0 function f_Ich(Ich)
  88. 0 Vth, g = hinv(Ich)
  89. 0 J = dhinv(Vth, g)
  90. 0 return abs(det(J)) * f_Vth_g(Vth, g)
  91. - end
  92. -
  93. 0 function f_Ich_pol(r1, ϕ1, r2, ϕ2)
  94. 1256640 Ich1 = pol2cart(r1, ϕ1)
  95. 1256640 Ich2 = pol2cart(r2, ϕ2)
  96. 1884960 Ich = @SMatrix[Ich1; Ich2]
  97. 58747920 return (r1 * r2) * f_Ich(Ich)
  98. - end
  99. -
  100. 0 function extremes(prob_interval, Finv_Vths, Finv_gs)
  101. 0 pts_Vths = [@SMatrix [Finv_Vths[1](p1); Finv_Vths[2](p2);;] for p1 = prob_interval for p2 = prob_interval]
  102. 0 pts_gs = [@SMatrix [Finv_gs[1](p1); Finv_gs[2](p2);;] for p1 = prob_interval for p2 = prob_interval]
  103. 0 [h(Vth, g) for Vth in pts_Vths for g in pts_gs]
  104. - end
  105. -
  106. 0 pts = extremes((0.0001, 0.9999), [Finv_Vth1, Finv_Vth2], [Finv_g1, Finv_g2])
  107. -
  108. 0 currents = reshape(collect(Iterators.flatten(pts)),2,length(pts))
  109. -
  110. 0 Ichs = mean(currents, dims=2)
  111. 0 Ich1 = Ichs[1]
  112. 0 Ich2 = Ichs[2]
  113. -
  114. 0 cart2pol(x) = (r=abs(x), ϕ=acos(real(x)/abs(x)) * sign(imag(x)))
  115. 0 pol2cart(r, ϕ) = r * (cos(ϕ) + sin(ϕ) * 1im)
  116. -
  117. 0 currents_pol = cart2pol.(currents)
  118. 0 extrema_r = extrema((x -> x.r) ∘ cart2pol, currents, dims=2)
  119. 0 extrema_ϕ = extrema((x -> x.ϕ) ∘ cart2pol, currents, dims=2)
  120. 0 int_r1, int_ϕ1 = extrema_r[1], extrema_ϕ[1]
  121. 0 int_r2, int_ϕ2 = extrema_r[2], extrema_ϕ[2]
  122. -
  123. 0 ich1_pol(r, ϕ) = f_Ich_pol(r, ϕ, cart2pol(Ich2)...) # fix ch2 current and get ch1
  124. 0 ich1_pol(r) = quadgk(ϕ -> ich1_pol(r, ϕ), int_ϕ1..., rtol=1e-3)[1]
  125. 0 ich1_pol_norm = quadgk(ich1_pol, int_r1..., rtol=1e-3)[1]
  126. - ich1(is) = ich1_pol(is) / ich1_pol_norm
  127. -
  128. 0 ich2_pol(r, ϕ) = f_Ich_pol(cart2pol(Ich1)..., r, ϕ) # fix ch1 current and get ch2
  129. 0 ich2_pol(r) = quadgk(ϕ -> ich2_pol(r, ϕ), int_ϕ2..., rtol=1e-3)[1]
  130. 0 ich2_pol_norm = quadgk(ich2_pol, int_r2..., rtol=1e-3)[1]
  131. - ich2(is) = ich2_pol(is) / ich2_pol_norm
  132. -
  133. - #return ich1, ich2
  134. - ## better integration
  135. -
  136. 0 f_Ich_pol_(x) = f_Ich_pol(x[1], x[2], x[3], x[4])
  137. 1884960 ich1_pol_(r1) = hcubature(x -> f_Ich_pol_([r1, x[1], x[2], x[3]]),
  138. - (int_ϕ1[1],int_r2[1],int_ϕ2[1]),
  139. - (int_ϕ1[2],int_r2[2],int_ϕ2[2]);
  140. - reltol=1e-3)[1]
  141. 0 ich2_pol_(r2) = hcubature(x -> f_Ich_pol_([x[2], x[1], r2, x[3]]),
  142. - (int_ϕ1[1],int_r1[1],int_ϕ2[1]),
  143. - (int_ϕ1[2],int_r1[2],int_ϕ2[2]);
  144. - reltol=1e-3)[1]
  145. -
  146. 0 return ich1_pol_, ich2_pol_
  147. - end
  148.  
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!