falkTX
/
Carla
mirror of https://github.com/falkTX/Carla
1
0
Fork 0
Code Releases 43 Activity
Audio plugin host https://kx.studio/carla
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1898 Commits
17 Branches
341MB
Tree: fdac1f6833
Branches Tags
${ item.name }
Create branch ${ searchTerm }
from 'fdac1f6833'
${ noResults }
Carla/source/modules/native-plugins/zamtube/wdf.cpp

658 lines
11KB
Raw Blame History 

  1. //wdf.cpp

  2. 

  3. #include <stdio.h>

  4. #include <inttypes.h>

  5. #include <cmath>

  6. #include "wdf.h"

  7. using std::abs;

  8. 

  9. WDF::WDF() {}

  10. 

  11. void WDF::setWD(T val) {

  12. 	WD = val;

  13. 	state = val;

  14. 	DUMP(printf("DOWN\tWDF\t%c\tWD=%f\tWU=%f\tV=%f\n",type,WD,WU,(WD+WU)/2.0));	

  15. }

  16. 

  17. void OnePort::setWD(T val) {

  18. 	WD = val;

  19. 	state = val;

  20. 	DUMP(printf("DOWN\tOneport\t%c\tWD=%f\tWU=%f\tV=%f\n",type,WD,WU,(WD+WU)/2.0));	

  21. }

  22. 

  23. T WDF::Voltage() {

  24. 	T Volts = (WU + WD) / 2.0;

  25. 	return Volts;

  26. }

  27. 

  28. T WDF::Current() {

  29. 	T Amps = (WU - WD) / (2.0*PortRes);

  30. 	return Amps;

  31. }

  32. 

  33. template <class Port1, class Port2>ser::ser(Port1 *l, Port2 *r) : Adaptor(THREEPORT) {

  34. 	left = l;

  35. 	right = r;

  36. 	PortRes = l->PortRes + r->PortRes;

  37. 	type = 'S';

  38. }

  39. 

  40. ser::ser(R* l, par* r) : Adaptor(THREEPORT) {

  41. 	left = l;

  42. 	right = r;

  43. 	PortRes = l->PortRes + r->PortRes;

  44. 	type = 'S';

  45. }

  46. 

  47. ser::ser(C* l, R* r) : Adaptor(THREEPORT) {

  48. 	left = l;

  49. 	right = r;

  50. 	PortRes = l->PortRes + r->PortRes;

  51. 	type = 'S';

  52. }

  53. 

  54. ser::ser(C* l, V* r) : Adaptor(THREEPORT) {

  55. 	left = l;

  56. 	right = r;

  57. 	PortRes = l->PortRes + r->PortRes;

  58. 	type = 'S';

  59. }

  60. 

  61. template <class Port>inv::inv(Port *l) : Adaptor(PASSTHROUGH) {

  62. 	left = l;

  63. 	PortRes = l->PortRes;

  64. 	type = 'I';

  65. }

  66. 

  67. inv::inv(ser *l) : Adaptor(PASSTHROUGH) {

  68. 	left = l;

  69. 	PortRes = l->PortRes;

  70. 	type = 'I';

  71. }

  72. 

  73. T ser::waveUp() {

  74. 	//Adaptor::WU = -left->waveUp() - right->waveUp();

  75. 	WDF::WU = -left->waveUp() - right->waveUp();

  76. 	DUMP(printf("UP\tser\tWU=%f\tWD=%f\tV=%f\n",WU,WD,(WD+WU)/2.0));

  77. 	return WU;

  78. }

  79. 

  80. template <class Port1, class Port2>par::par(Port1 *l, Port2 *r) : Adaptor(THREEPORT) {

  81. 	left = l;

  82. 	right = r;

  83. 	PortRes = 1.0 / (1.0 / l->PortRes + 1.0 / r->PortRes);

  84. 	type = 'P';

  85. }

  86. 

  87. par::par(inv* l, R* r) : Adaptor(THREEPORT) {

  88. 	left = l;

  89. 	right = r;

  90. 	PortRes = 1.0 / (1.0 / l->PortRes + 1.0 / r->PortRes);

  91. 	type = 'P';

  92. }

  93. 

  94. par::par(inv* l, V* r) : Adaptor(THREEPORT) {

  95. 	left = l;

  96. 	right = r;

  97. 	PortRes = 1.0 / (1.0 / l->PortRes + 1.0 / r->PortRes);

  98. 	type = 'P';

  99. }

  100. 

  101. par::par(C* l, R* r) : Adaptor(THREEPORT) {

  102. 	left = l;

  103. 	right = r;

  104. 	PortRes = 1.0 / (1.0 / l->PortRes + 1.0 / r->PortRes);

  105. 	type = 'P';

  106. }

  107. 

  108. T par::waveUp() {

  109. 	T G23 = 1.0 / left->PortRes + 1.0 / right->PortRes;

  110. 	WDF::WU = (1.0 / left->PortRes)/G23*left->waveUp() + (1.0 / right->PortRes)/G23*right->waveUp();

  111. 	DUMP(printf("UP\tpar\tWU=%f\tWD=%f\tV=%f\n",WU,WD,(WD+WU)/2.0));

  112. 	return WU;

  113. }

  114. 

  115. Adaptor::Adaptor(int flag) {

  116. 	WU = 0.0;

  117. 	WD = 0.0;

  118. 	switch (flag) {

  119. 		case ONEPORT:

  120. 			left = NULL;

  121. 			right = NULL;

  122. 			break;

  123. 		case PASSTHROUGH:

  124. 			right = NULL;

  125. 			break;

  126. 		default:

  127. 		case THREEPORT:

  128. 			break;

  129. 	}

  130. }

  131. 

  132. void ser::setWD(T waveparent) {

  133. 	Adaptor::setWD(waveparent);

  134. 	DUMP(printf("SER WP=%f\n",waveparent));

  135. 	left->setWD(left->WU-(2.0*left->PortRes/(PortRes+left->PortRes+right->PortRes))*(waveparent+left->WU+right->WU));

  136. 	right->setWD(right->WU-(2.0*right->PortRes/(PortRes+left->PortRes+right->PortRes))*(waveparent+left->WU+right->WU));

  137. }

  138. 

  139. void par::setWD(T waveparent) {

  140. 	Adaptor::setWD(waveparent);

  141. 	DUMP(printf("PAR WP=%f\n",waveparent));

  142. 	T p = 2.0*(waveparent/PortRes + left->WU/left->PortRes + right->WU/right->PortRes)/(1.0/PortRes + 1.0/left->PortRes + 1.0/right->PortRes);

  143. 

  144. 	left->setWD((p - left->WU));

  145. 	right->setWD((p - right->WU));

  146. }

  147. 

  148. T inv::waveUp() {

  149. 	///////////WD = -left->WD;

  150. 	WU = -left->waveUp(); 	//-

  151. 	DUMP(printf("UP\tinv\tWU=%f\tWD=%f\tV=%f\n",WU,WD,(WD+WU)/2.0));

  152. 	return WU;

  153. }

  154. 

  155. void inv::setWD(T waveparent) {

  156. 	WDF::setWD(waveparent);

  157. 	DUMP(printf("INV WP=%f\n",waveparent));

  158. 	//left->WD = -waveparent;		//-

  159. 	///////////left->WU = -WU;

  160. 	left->setWD(-waveparent);	//-

  161. 	

  162. }

  163. 

  164. R::R(T res) : Adaptor(ONEPORT) {

  165. 	PortRes = res;

  166. 	type = 'R';

  167. }

  168. 

  169. T R::waveUp() {

  170. 	WU = 0.0;

  171. 	DUMP(printf("UP\tR\tWU=%f\tWD=%f\tV=%f\n",WU, WD,(WD+WU)/2.0));

  172. 	return WU;

  173. }

  174. 

  175. C::C(T c, T fs) : Adaptor(ONEPORT) {

  176. 	PortRes = 1.0/(2.0*c*fs);

  177. 	state = 0.0;

  178. 	type = 'C';

  179. }

  180. 

  181. T C::waveUp() {

  182. 	WU = state;

  183. 	DUMP(printf("UP\tC\tWU=%f\tWD=%f\tV=%f\n",WU,WD,(WD+WU)/2.0));

  184. 	return WU;

  185. }

  186. 

  187. V::V(T ee, T r) : Adaptor(ONEPORT) {

  188. 	e = ee;

  189. 	PortRes = r;

  190. 	WD = 0.0;  //always?

  191. 	type = 'V';

  192. }

  193. 

  194. T V::waveUp() {

  195. 	T watts = 100.0;

  196. 	WU = 2.0*e - WD;

  197. 	if (Voltage()*Current() > watts) WU *= 0.999;//0.9955;

  198. 	DUMP(printf("UP\tV\tWU=%f\tWD=%f\tV=%f\n",WU, WD,(WD+WU)/2.0));

  199. 	return WU;

  200. }

  201. 

  202. inline T _exp(const T x)

  203. {

  204.     if(x < 10.0 && x > 10.0)

  205.         return 1.0 + x + x*x/2.0 + x*x*x/6.0 + x*x*x*x/24.0 + x*x*x*x*x/120.0

  206.             + x*x*x*x*x*x/720.0 + x*x*x*x*x*x*x/5040.0;

  207.     else

  208.         return exp(x);

  209. }

  210. 

  211. inline T _log(const T x)

  212. {

  213.     if(x > 10)

  214.         return log(x);

  215.     const T a=(x-1)/(x+1);

  216.     return 2.0*(a+a*a*a/3.0+a*a*a*a*a/5.0+a*a*a*a*a*a*a/7.0+a*a*a*a*a*a*a*a*a/9.0); 

  217. }

  218. 

  219. inline T _pow(const T a, const T b)

  220. {

  221.     return pow(a,b);

  222. }

  223. 

  224. T Triode::ffg(T VG) {

  225.         return (G.WD-G.PortRes*(gg*_pow(_log(1.0+_exp(cg*VG))/cg,e)+ig0)-VG);

  226. }

  227. 

  228. T Triode::fgdash(T VG) {

  229.         T a1 = exp(cg*VG);

  230.         T b1 = -e*pow(log(a1+1.0)/cg,e-1.0);

  231.         T c1 = a1/(a1+1.0)*gg*G.PortRes;

  232.         return (b1*c1);

  233. }

  234. 

  235. T Triode::ffp(T VP) { 

  236.     static bool prepared = false;

  237.     static double scale;

  238.     static double coeff[4];

  239.     if(!prepared) {

  240.         //go go series expansion

  241.         const double L2 = log(2.0);

  242. 

  243.         const double scale = pow(L2,gamma-2)/(8.0*pow(c,gamma));

  244.         coeff[0] = 8.0*L2*L2*scale;

  245.         coeff[1] = gamma*c*L2*4*scale;

  246.         coeff[2] = (c*c*gamma*gamma+L2*c*c*gamma-c*c*gamma)*scale;

  247.         coeff[3] = 0.0;

  248.         prepared = true;

  249.     }

  250. 

  251.     double A = VP/mu+vg;

  252.     return (P.WD+P.PortRes*((g*(coeff[0]+coeff[1]*A+coeff[2]*A*A))+(G.WD-vg)/G.PortRes)-VP);

  253. 

  254.     printf("%f\n", VP/mu+vg);

  255. 	return (P.WD+P.PortRes*((g*_pow(_log(1.0+_exp(c*(VP/mu+vg)))/c,gamma))+(G.WD-vg)/G.PortRes)-VP);

  256. }	//	    ^

  257. 

  258. T Triode::fpdash(T VP) {

  259.         T a1 = exp(c*(vg+VP/mu));

  260.         T b1 = a1/(mu*(a1+1.0));

  261.         T c1 = g*gamma*P.PortRes*pow(log(a1+1.0)/c,gamma-1.0);

  262.         return (c1*b1);

  263. }

  264. 

  265. T Triode::ffk() {

  266.         return (K.WD - K.PortRes*(g*pow(log(1.0+exp(c*(vp/mu+vg)))/c,gamma)));

  267. }

  268. /*

  269. T Triode::secantfg(T *i1, T *i2) {

  270.         T vgn = 0.0;

  271.         T init = *i1;

  272.         for (int i = 0; i<ITER; ++i) {

  273.                 vgn = *i1 - fg(*i1)*(*i1-*i2)/(fg(*i1)-fg(*i2));

  274.                 *i2 = *i1;

  275.                 *i1 = vgn;

  276.                 if ((fabs(fg(vgn))) < EPSILON) break;

  277.         }

  278.         if ((fabs(fg(vgn)) >= EPSILON)) {

  279.                 DUMP(fprintf(stderr,"Vg did not converge\n"));

  280. 		return init;

  281. 	}

  282.         return vgn; 

  283. }               

  284.         

  285. T Triode::newtonfg(T *i1) {

  286.         T init = *i1;

  287.         if (fabs(fg(*i1)) < EPSILON || fgdash(*i1)==0.0) return *i1;

  288.         T vgn = 0.0;

  289.         for (int i = 0; i<ITER; ++i) {

  290.                 vgn = *i1 - fg(*i1)/fgdash(*i1);

  291.                 *i1 = vgn;

  292.                 if (fabs(fg(vgn)) < EPSILON) break;

  293.         } 

  294.         if ((fabs(fg(vgn)) >= EPSILON)) {

  295. //                vgn = init;

  296.                 DUMP(fprintf(stderr,"Vg did not converge\n"));

  297.         }       

  298.         return vgn;

  299. }

  300. 

  301. T Triode::newtonfp(T *i1) {

  302.         T init = *i1;

  303.         if (fabs(fp(*i1)) < EPSILON || fpdash(*i1)==0.0) return *i1;

  304.         T vpn = 0.0;

  305.         for (int i = 0; i<ITER; ++i) {

  306.                 vpn = *i1 - fp(*i1)/fpdash(*i1);

  307.                 *i1 = vpn;

  308.                 if (fabs(fp(vpn)) < EPSILON) break;

  309.         }

  310.         if ((fabs(fp(vpn)) >= EPSILON)) {

  311. //                vpn = init;

  312.                 DUMP(fprintf(stderr,"Vp did not converge\n"));

  313.         }

  314.         return vpn;

  315. }

  316. 

  317. T Triode::secantfp(T *i1, T *i2) {

  318.         T vpn = 0.0;

  319.         for (int i = 0; i<ITER; ++i) {

  320.                 vpn = *i1 - fp(*i1)*(*i1-*i2)/(fp(*i1)-fp(*i2));

  321.                 *i2 = *i1;

  322.                 *i1 = vpn;

  323.                 if ((fabs(fp(vpn))) < EPSILON) break;

  324.         }

  325. 

  326.         if ((fabs(fp(vpn)) >= EPSILON))

  327.                 DUMP(fprintf(stderr,"Vp did not converge\n"));

  328.         return vpn;

  329. }

  330. */

  331. 

  332. //****************************************************************************80

  333. //	Purpose:

  334. //

  335. //		Brent's method root finder.

  336. //

  337. //	Licensing:

  338. //

  339. //		This code below is distributed under the GNU LGPL license.

  340. //

  341. //	Author:

  342. //

  343. //		Original FORTRAN77 version by Richard Brent.

  344. //		C++ version by John Burkardt.

  345. //		Adapted for zamvalve by Damien Zammit.

  346. 

  347. T Triode::r8_abs ( T x )

  348. {

  349. 	T value;

  350. 

  351. 	if ( 0.0 <= x )

  352. 	{

  353. 		value = x;

  354. 	}

  355. 	else

  356. 	{

  357. 		value = - x;

  358. 	}

  359. 	return value;

  360. }

  361. 

  362. Triode::Triode()

  363. {

  364. 	T r = 1.0;

  365. 

  366. 	while ( 1.0 < ( T ) ( 1.0 + r )	)

  367. 	{

  368. 		r = r / 2.0;

  369. 	}

  370. 

  371. 	r *= 2.0;

  372. 	r8_epsilon = r;

  373. }

  374. 

  375. T Triode::r8_max ( T x, T y )

  376. {

  377. 	T value;

  378. 

  379. 	if ( y < x )

  380. 	{

  381. 		value = x;

  382. 	}

  383. 	else

  384. 	{

  385. 		value = y;

  386. 	}

  387. 	return value;

  388. }

  389. 

  390. T Triode::r8_sign ( T x )

  391. {

  392. 	T value;

  393. 

  394. 	if ( x < 0.0 )

  395. 	{

  396. 		value = -1.0;

  397. 	}

  398. 	else

  399. 	{

  400. 		value = 1.0;

  401. 	}

  402. 	return value;

  403. }

  404. 

  405. 

  406. T Triode::zeroffp ( T a, T b, T t )

  407. {

  408. 	T c;

  409. 	T d;

  410. 	T e;

  411. 	T fa;

  412. 	T fb;

  413. 	T fc;

  414. 	T m;

  415. 	T macheps;

  416. 	T p;

  417. 	T q;

  418. 	T r;

  419. 	T s;

  420. 	T sa;

  421. 	T sb;

  422. 	T tol;

  423. //

  424. //	Make local copies of A and B.

  425. //

  426. 	sa = a;

  427. 	sb = b;

  428. 	fa = ffp ( sa );

  429. 	fb = ffp ( sb );

  430. 

  431. 	c = sa;

  432. 	fc = fa;

  433. 	e = sb - sa;

  434. 	d = e;

  435. 

  436. 	macheps = r8_epsilon;

  437. 

  438. 	for ( ; ; )

  439. 	{

  440. 		if ( abs ( fc ) < abs ( fb ) )

  441. 		{

  442. 			sa = sb;

  443. 			sb = c;

  444. 			c = sa;

  445. 			fa = fb;

  446. 			fb = fc;

  447. 			fc = fa;

  448. 		}

  449. 

  450. 		tol = 2.0 * macheps * abs ( sb ) + t;

  451. 		m = 0.5 * ( c - sb );

  452. 

  453. 		if ( abs ( m ) <= tol || fb == 0.0 )

  454. 		{

  455. 			break;

  456. 		}

  457. 

  458. 		if ( abs ( e ) < tol || abs ( fa ) <= abs ( fb ) )

  459. 		{

  460. 			e = m;

  461. 			d = e;

  462. 		}

  463. 		else

  464. 		{

  465. 			s = fb / fa;

  466. 

  467. 			if ( sa == c )

  468. 			{

  469. 				p = 2.0 * m * s;

  470. 				q = 1.0 - s;

  471. 			}

  472. 			else

  473. 			{

  474. 				q = fa / fc;

  475. 				r = fb / fc;

  476. 				p = s * ( 2.0 * m * q * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );

  477. 				q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );

  478. 			}

  479. 

  480. 			if ( 0.0 < p )

  481. 			{

  482. 				q = - q;

  483. 			}

  484. 			else

  485. 			{

  486. 				p = - p;

  487. 			}

  488. 

  489. 			s = e;

  490. 			e = d;

  491. 

  492. 			if ( 2.0 * p < 3.0 * m * q - abs ( tol * q ) &&

  493. 				p < abs ( 0.5 * s * q ) )

  494. 			{

  495. 				d = p / q;

  496. 			}

  497. 			else

  498. 			{

  499. 				e = m;

  500. 				d = e;

  501. 			}

  502. 		}

  503. 		sa = sb;

  504. 		fa = fb;

  505. 

  506. 		if ( tol < abs ( d ) )

  507. 		{

  508. 			sb = sb + d;

  509. 		}

  510. 		else if ( 0.0 < m )

  511. 		{

  512. 			sb = sb + tol;

  513. 		}

  514. 		else

  515. 		{

  516. 			sb = sb - tol;

  517. 		}

  518. 

  519. 		fb = ffp ( sb );

  520. 

  521. 		if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )

  522. 		{

  523. 			c = sa;

  524. 			fc = fa;

  525. 			e = sb - sa;

  526. 			d = e;

  527. 		}

  528. 	}

  529. 	return sb;

  530. }

  531. 

  532. T Triode::zeroffg ( T a, T b, T t )

  533. {

  534. 	T c;

  535. 	T d;

  536. 	T e;

  537. 	T fa;

  538. 	T fb;

  539. 	T fc;

  540. 	T m;

  541. 	T macheps;

  542. 	T p;

  543. 	T q;

  544. 	T r;

  545. 	T s;

  546. 	T sa;

  547. 	T sb;

  548. 	T tol;

  549. //

  550. //	Make local copies of A and B.

  551. //

  552. 	sa = a;

  553. 	sb = b;

  554. 	fa = ffg ( sa );

  555. 	fb = ffg ( sb );

  556. 

  557. 	c = sa;

  558. 	fc = fa;

  559. 	e = sb - sa;

  560. 	d = e;

  561. 

  562. 	macheps = r8_epsilon;

  563. 

  564. 	for ( ; ; )

  565. 	{

  566. 		if ( abs ( fc ) < abs ( fb ) )

  567. 		{

  568. 			sa = sb;

  569. 			sb = c;

  570. 			c = sa;

  571. 			fa = fb;

  572. 			fb = fc;

  573. 			fc = fa;

  574. 		}

  575. 

  576. 		tol = 2.0 * macheps * abs ( sb ) + t;

  577. 		m = 0.5 * ( c - sb );

  578. 

  579. 		if ( abs ( m ) <= tol || fb == 0.0 )

  580. 		{

  581. 			break;

  582. 		}

  583. 

  584. 		if ( abs ( e ) < tol || abs ( fa ) <= abs ( fb ) )

  585. 		{

  586. 			e = m;

  587. 			d = e;

  588. 		}

  589. 		else

  590. 		{

  591. 			s = fb / fa;

  592. 

  593. 			if ( sa == c )

  594. 			{

  595. 				p = 2.0 * m * s;

  596. 				q = 1.0 - s;

  597. 			}

  598. 			else

  599. 			{

  600. 				q = fa / fc;

  601. 				r = fb / fc;

  602. 				p = s * ( 2.0 * m * q * ( q - r ) - ( sb - sa ) * ( r - 1.0 ) );

  603. 				q = ( q - 1.0 ) * ( r - 1.0 ) * ( s - 1.0 );

  604. 			}

  605. 

  606. 			if ( 0.0 < p )

  607. 			{

  608. 				q = - q;

  609. 			}

  610. 			else

  611. 			{

  612. 				p = - p;

  613. 			}

  614. 

  615. 			s = e;

  616. 			e = d;

  617. 

  618. 			if ( 2.0 * p < 3.0 * m * q - abs ( tol * q ) &&

  619. 				p < abs ( 0.5 * s * q ) )

  620. 			{

  621. 				d = p / q;

  622. 			}

  623. 			else

  624. 			{

  625. 				e = m;

  626. 				d = e;

  627. 			}

  628. 		}

  629. 		sa = sb;

  630. 		fa = fb;

  631. 

  632. 		if ( tol < abs ( d ) )

  633. 		{

  634. 			sb = sb + d;

  635. 		}

  636. 		else if ( 0.0 < m )

  637. 		{

  638. 			sb = sb + tol;

  639. 		}

  640. 		else

  641. 		{

  642. 			sb = sb - tol;

  643. 		}

  644. 

  645. 		fb = ffg ( sb );

  646. 

  647. 		if ( ( 0.0 < fb && 0.0 < fc ) || ( fb <= 0.0 && fc <= 0.0 ) )

  648. 		{

  649. 			c = sa;

  650. 			fc = fa;

  651. 			e = sb - sa;

  652. 			d = e;

  653. 		}

  654. 	}

  655. 	return sb;

  656. }