+ /* some usefull expressions */
+ frac1=c2/(d2-h_cos2);
+ bracket1=1+c2d2-frac1;
+ if(f_c_ik==0.0) {
+ bracket2=0.0;
+ bracket2_n_1=0.0;
+ bracket2_n=0.0;
+ bracket3=1.0;
+ printf("Foo -> 0: ");
+ }
+ else {
+ bracket2=f_c_ik*bracket1;
+ bracket2_n_1=pow(bracket2,n-1.0);
+ bracket2_n=bracket2_n_1*bracket2;
+ bracket3=1.0+betan*bracket2_n;
+ printf("Foo -> 1: ");
+ }
+ bracket3_pow_1=pow(bracket3,(-1.0/(2.0*n))-1.0);
+printf("THETA: %.15f %.15f\n",cos_theta,theta*180/(2*M_PI));
+ bracket3_pow=bracket3_pow_1*bracket3;
+
+ /* now go on with calc of b_ij and derivation of b_ij */
+ b_ij=chi*bracket3_pow;
+
+ /* derivation of theta */
+ v3_scale(&force,&dist_ij,d_theta1);
+ v3_scale(&temp,&dist_ik,d_theta2);
+ v3_add(&force,&force,&temp);
+
+ /* part 1 of derivation of b_ij */
+ v3_scale(&force,&force,sin_theta*2*h_cos*f_c_ik*frac1);
+
+ /* part 2 of derivation of b_ij */
+ v3_scale(&temp,&dist_ik,df_c_ik*bracket1);
+
+ /* sum up and scale ... */
+ v3_add(&temp,&temp,&force);
+ scale=bracket2_n_1*n*betan*(1+betan*bracket3_pow_1)*chi*(1.0/(2.0*n));
+ v3_scale(&temp,&temp,scale);
+
+ /* now construct an energy and a force out of that */
+ v3_scale(&temp,&temp,f_a);
+ v3_scale(&force,&dist_ij,df_a*b_ij);
+ v3_add(&temp,&temp,&force);
+ v3_scale(&temp,&temp,f_c);
+ v3_scale(&force,&dist_ij,df_c*b_ij*f_a);
+ v3_add(&force,&force,&temp);
+