Jump to content


Check out our Community Blogs

Register and join over 40,000 other developers!


Recent Status Updates

View All Updates

Photo
- - - - -

Understand of code written in FOCAL

simulation focal

This topic has been archived. This means that you cannot reply to this topic.
2 replies to this topic

#1 cocowasabi

cocowasabi

    CC Lurker

  • Just Joined
  • Pip
  • 2 posts

Posted 07 June 2014 - 02:50 AM

Hey there! I come from a biology background so my programming experience is limited to basic Matlab. I'm trying to understand a simulation model of the aorta from 1971. The paper states that the simulation is written in standard FOCAL language, which seems so old, there is hardly any info on the web.

 

This is what I understand from the code so far:

- I understand the basic commands in the code (i.e, C, ASK, TYPE, SET)

- Therefore I understand how they have set up the basic equations for the simulation. 

 

What I have been trying to figure out:

- How does the integration command work?

- How are functions defined?

- The Control part (20.01) is a complete mystery to me.

02.01 C INTEGRATION CONTROL
02.10 TYPE !
02.12 ASK “TOTAL TIME”, TT
02.14 ASK “INTEGRATION INTERVAL”, DT
02.16 ASK ”(ITERATIONS BETWEEN) PRINTOUTS”, PI
02.18 ASK “(ITERATIONS BETWEEN) DISPLAYS”, SI
02.40 TYPE !!; DO 3; TYPE “SOLUTION”,!!
02.50 SET T=0; SET S0SI; SET P=PI
02.70 GOTO 4.01

04.01 C INITIAL CONDITIONS AND PARAMETERS
04.10 SET AV=80

05.01 C MODEL WITH INTEGRATORS FIRST
05.10 SET P=120*FSIN(453*T); C        DRIVING FUNCTION
05.20 SET VP=FUNC(P,0,0150,150)

06.10 SET AV=AV+NI*DT; C        INTEGRATOR
06.20 SET AP=FUNC(AV,0,0,50,20,100,70,200,220)
06.30 SET PF=FUNC(AP,20,0,50,1000,5000,200,15000)
06.40 SET PG=VP-AP
06.50 SET CV=FUNC(PG,0,0,0.00001,1000,200,1000)
06.60 SET VO=CV*PG
06.70 SET NI=VO-PF

20.01 C CONTROL
20.15 IF (SI - S) 20.20,20.20,20.21
20.20 SET S=1;00 22;GOTO 20.40
20.21 SET S=S+1
20.40 IF (PI-P) 20.42,20.42,20.44
20.42 SET P=1;00 23;GOTO 20.50
20.44 SET P=P+1
20.50 CONTINUE
20.54 IF (T-TT) 20.57,20.60,20.60
20.57 SET T=T+DT
20.58 GOTO 5.01
20.60 TYPE !!,”END OF SOLUTION”, !!
20.62 QUIT

Any help would be appreciated! Thanks



#2 0xDEADBEEF

0xDEADBEEF

    CC Devotee

  • Senior Member
  • PipPipPipPipPipPip
  • 790 posts

Posted 07 June 2014 - 11:52 AM

Found some references; 

http://homepage.cs.u...al/focal69.html

http://en.wikipedia.org/wiki/FOCAL-69

 

I think that there might be some library that the paper referneces that you've not posted - as the 'FUNC' command doesn't seem to be a standard command. Given what appears to be Random parameters to it.

 

The information in those links should be enough to convert this into pseduo code - let us know if you're still struggling with it.


Creating SEGFAULTs since 1995.


#3 cocowasabi

cocowasabi

    CC Lurker

  • Just Joined
  • Pip
  • 2 posts

Posted 09 June 2014 - 10:09 AM

Thank for you reply 0xDEADBEEF! I've been trying for a couple days now. I have looked at the references that you listed and tried to understand more.

 

So far I've managed to write out this "VP" driving function, which the code just defines as a truncated sine wave, it looks spot on, like in the paper:

%%ventricular pressure
t=0:0.001:0.02;
VP= zeros(size(t));
for i=1:length(t)
    if 0.007<=t(i) && t(i)<0.014
        VP(i) = 0;
    else VP(i) = 120*sin((6.28/0.0139)*t(i));
    end
end

plot(t,VP);
grid on;
hold on;

The paper just listed FUNC as a "graphical function generator," which I didn't understand if MATLAB had something similar, which just defines some arbitrary function with the given parameters.

 

I found another reference which defines the 'AP' and 'PF' FUNC's with empirical equations. I used the parameters listed above to find the missing coefficients and now I have those equations.

AP=(1/30).*AV + (1/125).*(AV.^2) + (-1/75000).*(AV.^3);
PF=(500/3) +(-55/2).*AP + (121/120).*(AP.^2) + (-1/400).*(AP.^3);

Other than that, I've just typed in the other equations into MATLAB.

%% initial conditions
AV=80;
CV=1000;

%% equations
NI=VO-PF;
AV=AV+NI*t; %integral of NI

% VO=CV*(VP-AP); but only when (VP-AP)>0 else VF=0

VO= zeros(size(t));
for j=1:length(t)
    if 0.007<=t(j) && t(j)<0.014
        VO(j) = 0;
    else VO(j) = ((120*sin(453*t(j)))-AP).*CV;
    end
end

I don't really know what I was doing with the second 'for' loop. I just kinda copied the first because it worked. I know that for VO, I want to only consider the positive solutions, so I thought that copying the first 'for' loop was how to do it, but I'm not sure.

 

I can't see any other graphical results to compare with the paper because all the equations seem to depend on each other, and I'm going around in an endless loop!

In the end I trying to get a graphical output for VP(check!), VO, and AP with respect to time.

 

Specifically, what I'm trying to figure out is:

- Why does the original code say AV=80 as the initial parameter? How can AV=AV+NI*t?? Shouldn't it be AV_i (AV initial)=80? But not AV itself? Does this change anything?

- How do I decipher these equations that all depend on each other? Is this an ode, should I be using some kind of ode solver, some integration function?

 

Please, excuse my ignorance, i'm trying, but utterly confused.

Thanks for helping me,

cocowasabi






Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download