physics engine - phase order and other general inf

I want to build my own 2D mini-physics engine, that will include(for now) rigid bodys, and constraints (joints, contacts, springs...). And I've tried to figure out what is the right order of the phases, when I start a timestep, when the the general phases are: Broadphase, Narrow phase(Collision Detection, and Contact generation), Resolution (Constraints Solver), and the Integration - hope you can tell me the right order.

I allso have general questions about each phase:

Narrow Phase - When I've found a collision do I need to seperate the bodies after I've found the collision or just apply an impulse to their velocities in the resolution phase?

and if I'm using CCD(Continuous Collision Detection) for some bodies, what happand if I've found a collision (the earliest for object A) of a fast moving object A with other object B, and than I found that object B will collide with other object C before object he will collide with object A, do I need to go back to object A and looking for other collisions after the first collision I've found to him?

and if I'm using Contact Solver in my engine (in the constraints solver, I'll generat A Contact constraints for each collision, and I'll solve then all in the resolution phase, or I'll solve each contact right after I've found them in the narrow phase?

Resolution - Using a constaints Solver for the resolution phase is ok? and what if by solving one constraint I'll solve other constraints or make more constraints? (Do I need to check after each constraints solving if it solves other constraints or make more?

Intefration - this is the part where I take all the information I gathered (impulses, forces..) and intefrate them to velocity and than to position with the integration method I chose?

A symplectic euler integration will be enough for my physics engine?

and I alsso saw in many physics engines, like box2D, they use iterations and let me chose the frequncy like here (notice I can change the iteration count (10) and the frequency (60.0 Hz)):

What does this variables means? the iterations is how many times I recall the physics update(all the phases above) in each frame, or it just recall one phase like the narrow phase or somthing like that? and the frequency variable let me chose how many frames there will be in one second? correct me if i'm worng please.

all the qeustions above get me stack from learning because its basics things, and i didnt found a place that explain this things in a clean and strait forward

so thanks for anyone who will read all of that, and especially for anyone who will also help me with all of my questions :)

the order of your phases are not so important. Yes they do mess each other but any dis-function can be corrected by adjusting constants ...

getting started

you should start with object representation and visualization prior to any simulation code. You need to be able to visualize your stuff so you can see if is all OK or not
object classes

so create class(es) for all object types you want to support and add few basic interface functions (virtual base class is a good idea for future implementation) like:
- load(file,ini,stream...),save(file,ini,stream...)
- draw(screen or render context),update(dt),bool colide(object),etc...
in time you will see exactly what you need to add ... They will be empty for now (instead the draw which you need to code...). You can also prepare few common physics variables (also as virtual) like position, speed, temperature,...
base class

it should be the encapsulation of your engine. It hold the lists of all the objects/stuff. Main interface with GUI/App. For more insight of what I mean by now look here Drag&Drop editor in C++ you can use it as start point just add the physics iteration/update to it. Do not forget to implement the object interface into this main class like:
- load(filename) ... which load all the objects
- draw(...) ... which draws all the objects ...
simulation

start with this only if the above is already working. You need to take into account the accuracy and response time you want to achieve. The more fast object/processes are the higher simulation loop frequency you need. For mine simulations I usually use 1-100 ms per iteration. For very fast simulation you can do N iterations per one timer call. That is the meaning of the count and frequency. Where frequency is the timer call speed your iteration loop is called and count (N) is the time step dt division so if N=10 and f=60Hz it means you iterate 10 times every 1/60 seconds so it is the same as N=1 f=600Hz but at least in Windows timer resolution is 1 ms and not very precise ether. So frequencies above 100 Hz are unreliable. if you want to be even more precise you can measure the time by yourself via PerformanceCounter or RDTSC or any precise enough time API you have at disposal.

you can use D'Lambert principle (simple integration) for the simulation of motions the rest compute with known equations (do not know all of what you want to simulate). For example look here simple mass points gravity simulation in 3D (C++)

if your code is well written you can parallelize the iteration loop but you have to take in mind that it can create few troubles like double collision reaction,etc
collisions

there are tons of stuff about this here on SO/SE so just search for it. I think you should look here: simulation of particle collisions to get the idea of multiple collisions handling.

To get the idea here is example of how one of mine simulations looks like

it uses all the stuff I wrote here about and in addition it uses a special class for bonds like springs or joints. This is how it looks like:

struct _bond
    {
    physics_point *pnt0,*pnt1;
    double l0,l1;
    int _beg0,_end0;
    int _beg1,_end1;
    List<int> depend0,depend1;
    int _computed;
    _bond()     {}
    _bond(_bond& a) { *this=a; }
    ~_bond()    {}
    _bond* operator = (const _bond *a) { *this=*a; return this; }
    //_bond* operator = (const _bond &a) { ...copy... return this; }
    };
List<_bond> bnds;

I could cover whole book with the stuff behind this but I am too lazy for that and also this site is not the right place for this so just in a hurry the important stuff is that each bond has it start and end object pointer (pnt0,pnt1) and during computations are filled the dependency lists (depend0,depend1).

depend0 for all consequential bonds from pnt0 side
depend1 for all consequential bonds from pnt1 side.

Then during each iteration is recursively updated all bonds to match the bond/collision conditions only when they all do then the positions are updated at once for all of them.