RK.h (1409B)
1 #ifndef RK_H 2 #define RK_H 3 4 bool RK_step(double state[]){ 5 //state={q1,q2,p1,p2} 6 //if open orbit step won't compute 7 if (state[0]>10 || state[1]>10 || state[0]<-10 || state[1]<-10){ 8 return false; 9 } 10 double k1[4]={0,0,0,0}; 11 double k2[4]={0,0,0,0}; 12 double l1[4]={0,0,0,0}; 13 double l2[4]={0,0,0,0}; 14 k1[0]=dt*state[2]; 15 k2[0]=dt*state[3]; 16 l1[0]=dt*(-state[0]-2*state[0]*state[1]); 17 l2[0]=dt*(-state[1]-state[0]*state[0]+state[1]*state[1]); 18 k1[1]=dt*(state[2]+l1[0]/2); 19 k2[1]=dt*(state[3]+l2[0]/2); 20 l1[1]=dt*(-(state[0]+k1[0]/2)-2*(state[0]+k1[0]/2)*(state[1]+k2[0]/2)); 21 l2[1]=dt*(-(state[1]+k2[0]/2)-(state[0]+k1[0]/2)*(state[0]+k1[0]/2)+(state[1]+k2[0]/2)*(state[1]+k2[0]/2)); 22 k1[2]=dt*(state[2]+l1[1]/2); 23 k2[2]=dt*(state[3]+l2[1]/2); 24 l1[2]=dt*(-(state[0]+k1[1]/2)-2*(state[0]+k1[1]/2)*(state[1]+k2[1]/2)); 25 l2[2]=dt*(-(state[1]+k2[1]/2)-(state[0]+k1[1]/2)*(state[0]+k1[1]/2)+(state[1]+k2[1]/2)*(state[1]+k2[1]/2)); 26 k1[3]=dt*(state[2]+l1[2]); 27 k2[3]=dt*(state[3]+l2[2]); 28 l1[3]=dt*(-(state[0]+k1[2])-2*(state[0]+k1[2])*(state[1]+k2[2])); 29 l2[3]=dt*(-(state[1]+k2[2])-(state[0]+k1[2])*(state[0]+k1[2])+(state[1]+k2[2])*(state[1]+k2[2])); 30 state[0]+=(k1[0]+2*k1[1]+2*k1[2]+k1[3])/6; 31 state[1]+=(k2[0]+2*k2[1]+2*k2[2]+k2[3])/6; 32 state[2]+=(l1[0]+2*l1[1]+2*l1[2]+l1[3])/6; 33 state[3]+=(l2[0]+2*l2[1]+2*l2[2]+l2[3])/6; 34 return true; 35 } 36 37 #endif