B.1 Simulationsprogramme zur Kontrastinversion

An dieser Stelle ist das im Rahmen dieser Diplomarbeit erstellte Programm zur Simulation der Kontrastinversion mit und ohne Abschattungseffekt aufgelistet. Es wird das hexagonale harmonische Potential, das in Anhang D abgebildet ist, verwendet. Bei den Programmen zu den anderen Simulationen, insbesondere dem Programm zur Verdeutlichung der abgeschatteten Bereiche, handelt es sich um Modifikationen auf der Basis des folgenden Programms, die aus Platzgründen nicht im einzelnen abgedruckt werden.

/* sim.h

Author: Peter Münster, 1996

Kopf für alle Module

*/

#if !defined( __SIM__ )
#define __SIM__

// #include <curses.h>
#include <string.h>
#include <unistd.h>
#include <stdio.h>
#include "/home/peter/c/mystd.h"
#include <math.h>

/*
#pragma warn -pia
#pragma warn -rvl
*/

typedef unsigned unsgnd;

typedef struct proc
{
pid_t pid;
char name[DAT_NAME_MAX];
}
proc_struct;

typedef struct info
{
char doc[DAT_NAME_MAX]; // Info-String
// char dir[DAT_NAME_MAX]; // current working directory
doub phi; // "Barrierenhöhe"
doub cut; // Abbruchkriterium
doub ant; /* Anteil der nichtabgeschatteten Punkte im
   letzten Kreisring */
int rmax; // maximaler Radius für Abschattung in Pixeln
doub scal; // Skalierungsfaktor = Pixel pro Ångstrøm
int dim; // Breite = Höhe in Pixeln
doub c[2]; // Gitterkonstante und Korrugation in Ångstrøm
doub d0; // Abstand in Ångstrøm
doub dd0; // Schrittweite
doub d0max; // maximaler Abstand in Ångstrøm
int flag; // Schalter: 0=normal, 1=genau0, 2=genau1
}
info_struct;

#define MAX(x,y) ((x)>(y)?(x):(y))

// 3 nützliche Makros falls curses benutzt wird
#define pf1(x,y) printf(x,y);fflush(NULL)
#define pf2(x,y,z) printf(x,y,z);fflush(NULL)
#define mywrite(x) write(1, x, strlen(x))

// lese Zeichen mit oder ohne Echo
int getche(void); // in simhilf.c
int getchne(void); // in simhilf.c

// von Max
void schreibedim(int dim, FILE *ausdp); // in simhilf.c

// gibt max bzw. min aus Array zurück
doub maximum(doub *val, int j); // in simhilf.c
doub minimum(doub *val, int j); // in simhilf.c
// weil in math.h kein round ist:
doub round(doub x); // in simhilf.c

/*liest globale Info-Struktur aus datei,
  oder schreibt sie in datei */
int ladeinf(char *datei, char *mode,
size_t (*func)(void *ptr,
   size_t size,
   size_t nmemb,
   FILE *stream)); // in siminfo.c
// zeigt Info-Struktur an
int zeigeinfo(char *datei); // in siminfo.c

// rechnet und erzeugt die Ausgabedateien
void rechne(char *kontrastfile); // in simcalc.c

// kleine Hilfsfunktion für die Eingabe
void my_doub_scan(char *s, doub *d_ptr); // in siminput.c

// Eingabe der Info-Struktur per Hand
int eingabe(void); // in siminput.c

// noch zu erledigen
int ausgabe(char *datei); // in simoutput.c

// durchrechnen von Serien
void serie(proc_struct **pr_ptr, int *cnt); // in simserie.c

// Hilfsfunktion für serie
void serienfunc(char *str); // in simserie.c

#endif

Das anschließende C-Modul enthält die Hauptfunktion main, in der bereits zusätzlich einige Menüpunkte wie etwa bestimmte Ausgaberoutinen vorgesehen sind, die jedoch nicht mehr implementiert wurden.

/* simmain.c

Author: Peter Münster, 1996

Hauptmodul für die Simmulation des Tunnelstromes beim STM

*/

#include "sim.h"
#include <signal.h>

#if DEBUG
#define FORK_UND_RECHNE(xxx) /* nur Textersatz */ \
fork_ret=fork(); /* wird für wait gespeichert */\
switch(fork_ret){ \
case -1: \
FPF("fork-Aufruf in main %s", "ohne Erfolg!\n");\
break; \
case 0: \
exit(0); /* ohne Interesse für’s Debuggen */\
forkcount++; /* Prozeßinfos abspeichern: */ \
proc_ptr=(proc_struct *)realloc(proc_ptr,\
sizeof(proc_struct)*forkcount); \
proc_ptr[forkcount-1].pid=fork_ret; \
strncpy(proc_ptr[forkcount-1].name,\
xxx, DAT_NAME_MAX); \
proc_ptr[forkcount-1].name[DAT_NAME_MAX-1]=’\0’;\
default: \
rechne(""); /* rechne = Vaterprozeß */    \
exit(0); \
}
#else

#define FORK_UND_RECHNE(xxx) /* nur Textersatz */ \
fork_ret=fork(); /* wird für wait gespeichert */\
switch(fork_ret){ \
case -1: \
FPF("fork-Aufruf in main %s", "ohne Erfolg!\n");\
break; \
case 0: \
rechne(""); \
exit(0); \
default: \
forkcount++; /* Prozeßinfos abspeichern: */ \
proc_ptr=(proc_struct *)realloc(proc_ptr,\
sizeof(proc_struct)*forkcount); \
proc_ptr[forkcount-1].pid=fork_ret; \
strncpy(proc_ptr[forkcount-1].name,\
xxx, DAT_NAME_MAX); \
proc_ptr[forkcount-1].name[DAT_NAME_MAX-1]=’\0’;\
}

#endif

char progname[DAT_NAME_MAX];

extern info_struct inf;

void finish(int sig){
char str[DAT_NAME_MAX];
sprintf(str, ".simrc.%2.0lf", inf.cut);
ladeinf(str, "w", fwrite);
sleep(3);
#if ABSCHATTUNG
sprintf(str,
"echo nice %s .simrc.%2.0lf|at now+2hours 2>/dev/null",
"/home/peter/bin/simma", inf.cut);
#else
sprintf(str,
"echo nice %s .simrc.%2.0lf|at now+2hours 2>/dev/null",
"/home/peter/bin/simoa", inf.cut);
#endif
system(str);

#if DEBUG
FPF("finish: erhielt Signal Nummer %d!\n", sig);
#endif
exit(0);
}

int main(int argc, char **argv) {
int i, c=1, forkcount=0;
pid_t fork_ret;
char str[DAT_NAME_MAX];
FILE *fp; /* in proc_struct werden Prozeßinfos
   gespeichert */
proc_struct *proc_ptr=(proc_struct *)
malloc(sizeof(proc_struct));

signal(SIGTERM, finish); // für die Automatik

strncpy(progname, argv[0], DAT_NAME_MAX-1); /* progname ist
   global!
   (für Fehlermeldungen) */

ladeinf(".simrc", "r", fread); /* inf-Struktur laden
   (enthält
   Probeninformationen) */

if(argc==2){ // falls ein Parameter, gehe auf Automatik!
ladeinf(argv[1], "r", fread);
for(; inf.d0<inf.d0max; inf.d0+=inf.dd0)
rechne("");
ladeinf(argv[1], "w", fwrite);
return 0;
}

#if CURSES
initscr(); /* initialize the curses library */
nocbreak(); /* take input chars one at a time,
   no wait for \n */
echo();
#endif

if(argc>2){ // falls mehrere Parameter, kein In- und Output
while(--argc){
if(ladeinf(argv[argc], "r", fread)) /* lädt die
inf-Struktur */
continue; // Fehler!
FORK_UND_RECHNE(argv[argc]);
}
return 0; /* Programmende, auf Sohnprozesse
   wird NICHT gewartet */
}

// system("clear"); // Clear-Home
pf1("\t\t\tSTM-Simulation\n%s",
"\t\t\t--------------\n\n"); // pf1: in sim.h
mywrite(
"Dateien werden in im Arbeitsverzeichnis gesucht!\n\n");

ladeinf(".simrc", "r", fread); /* inf-Struktur laden
   (enthält Probeninformationen) */
while (c) {
switch (c) {
case 1:
while(c){
mywrite("\n Eingabe\n");
mywrite(  "    ---------\n\n");
mywrite(" 0.) Keine\n");
mywrite(" 1.) Per Hand\n");
mywrite(" 2.) Aus Datei laden\n");
mywrite(" 3.) Serie\n\n ?  ");
scanf("%d", &c);
getchar(); /* das Returnzeichen bleibt
   sonst im Puffer */
switch(c){
case 1:
if(!eingabe()){ /* lädt die inf-Struktur
   neu (0 bei Erfolg) */
pf1(
"\nAls %s.inf im Verzeichnis daten speichern? (j/n)",
inf.doc);
scanf("%c", str); // "Mißbrauch" von str
if(*str != ’n’){
ladeinf(inf.doc, "w", fwrite);
mywrite("\nRechnen? (j/n)");
scanf("%s", str); /* "Mißbrauch"
   von str */
if(*str == ’n’)
break;
}
FORK_UND_RECHNE(inf.doc); /* wenn nicht
gespeichert wird, wird auf jeden Fall gerechnet */
}
break;
case 2:
mywrite("Aus welcher? (ohne .inf!)\n");
system("ls *.inf");
mywrite("\n? ");
scanf("%s", str);
str[DAT_NAME_MAX-1]=’\0’;// zur Sicherheit
if(!ladeinf(str, "r", fread))
FORK_UND_RECHNE(inf.doc);
break;
case 3:
if(!eingabe())
serie(&proc_ptr, &forkcount);
// Parameter dienen
break; // nur dem Aktualisieren
} // der Prozeßinformationen
}
case 2:
while(c){
mywrite("    Doku erzeugen\n");
mywrite("   ---------------\n\n");
mywrite(" 0.) Keine\n");
mywrite(" 1.) Die aktuelle\n");
mywrite(" 2.) Andere\n");
scanf("%d", &c);
switch(c){
case 1:
//makedocu(); noch zu erledigen
break;
case 2:
system("ls *.inf");
mywrite("\n ?  ");
scanf("%s", str);
str[DAT_NAME_MAX-1]=’\0’;
/*if(! ladeinf(str, "r", fread))
makedocu();*/
break;
}
}
case 3:
while(c){
mywrite("   Ausgabe\n");
mywrite(" ---------\n\n");
mywrite(" 0.) Keine\n");
mywrite(" 1.) Eine\n");
mywrite(" 2.) Mehrere\n");
scanf("%d", &c);
switch(c){
case 1:
system("ls x*");
mywrite("\n? ");
scanf("%s", str);
str[DAT_NAME_MAX-1]=’\0’;
//ausgabe(str); noch zu erledigen
break;
case 2:
#if 0
do{ noch zu erledigen
system("ls x*");
mywrite("\n(Return für Ende) ? ");
scanf("%s", str);
str[DAT_NAME_MAX-1]=’\0’;
}while(! ausgabe(str));
#endif
break;
}
}
break;
default:
mywrite("\n Und nochmal!\n");
}
mywrite("\n\n Menu:\n");
mywrite(" 0.) Ende\n");
mywrite(" 1.) Eingabe\n");
mywrite(" 2.) Doc-Datei erzeugen\n");
mywrite(" 3.) Ausgabe\n\n ? ");
scanf("%d", &c);
}
// endwin(); // falls curses benutzt wird
getchar(); // wegen letztem scanf
ladeinf(".simrc", "w", fwrite);
pf1(
"\nSoll auf die %d noch laufenden Prozesse gewartet werden? ",
forkcount);
scanf("%c", &c);
#if ! HP
if(c!=’n’) // OK, es wird gewartet
while(forkcount--){
pid_t pid=wait(&fork_ret);
for(i=0; proc_ptr[i].pid != pid; i++);
printf("\n%s: wait-Status: | %x | %x |\n",
   proc_ptr[i].name,
   (fork_ret>>8) & 0xff, fork_ret & 0xff);
}
#endif
return 0;
}

Die folgenden Dateien enthalten jeweils Funktionen für die Eingabe, das Schreiben und Lesen von Parametern, die eigentliche Simulationsrechnung, die Unterstützung von Simulationsserien sowie weitere nützliche Funktionen.

Zu dem Algorithmus für die Abschattung ist zu erwähnen, daß dazu Linien erzeugt werden, die daraufhin geprüft werden, ob sie komplett oberhalb der Oberfläche liegen oder nicht. Um Rechenzeit zu sparen —was wichtig ist, da diese Überprüfung innerhalb der innersten Schleife stattfindet—, wurden die bereits berechneten Punkte zur Linienerzeugung verwendet, was dazu führt, daß diese nicht mehr gerade, sondern treppenförmig wird, falls sie nicht auf der x- oder y-Achse liegt. Dies hatte zur Folge, daß die Simulationsbilder keine hexagonale Symmetrie mehr aufwiesen, denn diagonale Linien wurden schlechter aufgelöst als waagerechte oder senkrechte. Um diesem Fehler zu begegnen, durften Linien nicht mehr gerastert werden, sondern mußten nun für jeden Punkt analytisch ausgewertet werden, was zu dem Makro „SURFUNC(x,y)“ führte.

/* siminput.c

Author: Peter Münster, 1996

Eingabemodul für simmain.c

*/

#include "sim.h"

extern info_struct inf;

char progname[DAT_NAME_MAX];

// schreibt nach dp neuen Wert, oder läßt Default-Wert

void my_doub_scan(char *s, doub *dp){
char str[DAT_NAME_MAX+1];
pf1(s, *dp);
if(*((fgets(str, DAT_NAME_MAX, stdin))+1)) // \n ist min-
// destens das
*dp=(doub)strtod(str, NULL); // erste Zeichen (-> ..+1)
}

// Hauptfunktion für die Eingabe (noch ohne Fehlerbehandlung)

int eingabe(){
char str[DAT_NAME_MAX+1];

my_doub_scan(
"\nBitte die Barrierenhöhe eingeben!\n(default=%lf): ",
&inf.phi);
my_doub_scan("\nWann ist Abbruch?\n(default=%lf): ",
&inf.cut);
printf("\nBitte Modus eingeben!\n");
pf1("(default=%d):\n0=normal, 1=genau0, 2=genau1: ",
inf.flag);
if(*((fgets(str, DAT_NAME_MAX, stdin))+1)) // \n ist
// mindestens das
inf.flag=atoi(str); // erste Zeichen (-> ..+1)

/* my_doub_scan(
"\nBitte die Gewichtung eingeben!\n(default=%lf): ",
&inf.g); */
my_doub_scan(
"\nSkalierungsfaktor (Pixel/Å) eingeben!\n(default=%lf): "
,&inf.scal);
/* pf1("\nBitte rmax eingeben!\n(default=%d): ", inf.rmax);
if(*((fgets(str, DAT_NAME_MAX, stdin))+1)) // \n ist
// mindestens das
inf.rmax=atoi(str); // erste Zeichen (-> ..+1)
pf1(
   "\nBitte Breite (=Höhe) in Pixel eingeben!\n(default=%d): ",
   inf.dim);
   if(*((fgets(str, DAT_NAME_MAX, stdin))+1)) // \n ist
// mindestens das
inf.dim=atoi(str); // erste Zeichen (-> ..+1)
*/
my_doub_scan(
"\nBitte Gitterkonstante eingeben!\n(default=%lf): ",
&inf.c[0]);
my_doub_scan(
"\nBitte Korrugation eingeben!\n(default=%lf): ",
&inf.c[1]);
my_doub_scan("\nBitte Abstand eingeben!\n[%lf]: ", &inf.d0);

sprintf(inf.doc,
"%1.1lf#%1.6lf#%1.1lf#%1.1lf#%1.1lf#%1.1lf",
inf.phi,inf.cut,inf.scal,inf.c[0],
inf.c[1],inf.d0);

#if 0
pf2(
"\nBitte Doku bis zu %d Zeichen eingeben!\n(default=%s): ",
DAT_NAME_MAX, inf.doc);
if(*((fgets(str, DAT_NAME_MAX, stdin))+1)){ // wegen Return
strcpy(inf.doc, str);
inf.doc[strlen(str)-1]=’\0’;
}
#endif
return 0;
}

/* siminfo.c

Author: Peter Münster, 1996

Hilfsmodul für simmain.c

*/

#include "sim.h"

info_struct inf;

extern char progname[];

// lädt die globale inf-Struktur

int ladeinf(char datei[DAT_NAME_MAX], char mode[3],
size_t (*func)(void *ptr, size_t size, size_t nmemb,
FILE *stream)){
FILE *fp;

if(!(fp=fopen(datei, mode)))
return 1;
if(!func((void *)(&inf.phi), SOD, 1, fp))
return 1;
if(!func((void *)(&inf.cut), SOD, 1, fp))
return 1;
if(!func((void *)(&inf.scal), SOD, 1, fp))
return 1;
if(!func((void *)(&inf.c), SOD, 2, fp))
return 1;
if(!func((void *)(&inf.d0), SOD, 1, fp))
return 1;
if(!func((void *)(&inf.d0max), SOD, 1, fp))
return 1;
if(!func((void *)(&inf.dd0), SOD, 1, fp))
return 1;
if(!func((void *)(&inf.flag), SOI, 1, fp))
return 1;
fclose(fp);
return 0;
}

/*zeigt Indodatei in lesbarem Format und
füllt die globale Struktur inf */

int zeigeinfo(char datei[DAT_NAME_MAX]){
pf1("\nInhalt von %s.inf:\n\n", datei);
if(ladeinf(datei, "r", fread)){
FPF("\nzeigeinfo: aus %s.inf konnte nicht gelesen werden!",
datei);
return 1;
}
pf1("Barrierenhöhe: %lf\n", inf.phi);
pf1("Skalierungsfaktor (Pixel/Å): %lf\n", inf.scal);
// pf1("Breite und Höhe in Pixeln: %d\n", inf.dim);
pf1("Gitterkonstante: %lf\n", inf.c[0]);
pf1("Korrugation: %lf\n", inf.c[1]);
// pf1("Atomradius: %lf\n", inf.c[4]);
pf1("Abstand: %lf\n\n", inf.d0);
return 0;
}

/* simcalc.c

Author: Peter Münster, 1996

Rechenmodul für sim

*/

#include "sim.h"

extern char progname[];
extern info_struct inf;

static doub abstand; // vom Tal gemessen
static doub git; // in Pixel
static doub kor; // in Pixel
static doub *probfeld;
static short *flagfeld;
static doub sqrphi_scal; // für die exp-Funktion
static doub sq3; // für das Makro surfunc
static doub hilf; // für das Makro surfunc

// Erzeugen der Oberfläche
void surface(doub *probfeld){
register int i, j;
doub w0[2], w1[2], w2[2];

// die drei w#’s sind die reziproken Gittervektoren:
w0[0]=0 *2*M_PI/git;
w0[1]=1 *2*M_PI/git;
w1[0]=-sqrt(3)/2 *2*M_PI/git;
w1[1]=-0.5 *2*M_PI/git;
w2[0]=sqrt(3)/2 *2*M_PI/git;
w2[1]=-0.5 *2*M_PI/git;

// Hexagonale Cosinusfunktion, siehe Chen: Seite 132
for(j = 0; j < inf.dim; ++j)
for(i = 0; i < inf.dim; ++i){
probfeld[inf.dim*j + i] = 1.5;
probfeld[inf.dim*j + i] += cos(i*w0[0]+j*w0[1]);
probfeld[inf.dim*j + i] += cos(i*w1[0]+j*w1[1]);
probfeld[inf.dim*j + i] += cos(i*w2[0]+j*w2[1]);
probfeld[inf.dim*j + i] *= kor*2.0/9.0;
}
}

/*
// diese Funktion berechnet eine Linie zwischen (x,y) und
// (x+dy,y+dy)
void linie(int max, int dx, int dy, int x, int y,
  point_struct *line){
register int i;
for(i=0; i<max; i++){
line[i].x=x+rint((doub)dx*(doub)i/(doub)max);
line[i].y=y+rint((doub)dy*(doub)i/(doub)max);
}
}
*/

// surfunc als Makro aus Geschwindigkeitsgründen

#define SURFUNC(x,y)  ((1.5+cos((y)*2*hilf)+cos((x)*sq3*  \
hilf+(y)*hilf)+cos((x)*sq3*hilf-(y)*  \
   hilf))*kor*2/9)
// nützliche Abkürzungen:
#define LINIE_XY(xy,dxy,m,i) ((xy)+(doub)(dxy)* \
(doub)(i)/(doub)(m))
#define LINIE_Z(x,y,dx,dy,m,i) (SURFUNC(LINIE_XY(x,dx,m,i), \
LINIE_XY(y,dy,m,i)))

// getstrom berechnet den Strom von (i,j) nach (x,y)
doub getstrom(int x, int y, int i, int j){
doub d, wert=probfeld[j*inf.dim+i]; // wert= Funktionswert
doub index; // am Aufpunkt
int dx=i-x;
int dy=j-y;
doub dist=sqrt(dx*dx+dy*dy);

#if ABSCHATTUNG
index=(abstand-1.01*kor)*dist/(abstand-wert);
/* 1.01 soll nur einen Sicherheitsabstand gewährleisten */
while(index < dist && abstand-index*(abstand-wert)/dist >
  LINIE_Z(x,y,dx,dy,dist,index))
index++;
if(dist<=index){ // keine Abschattung
#endif
flagfeld[j*inf.dim+i]=2; // 2: schon berechnet

if(inf.flag)
d=sqrt(dx*dx+dy*dy+(abstand-wert)*(abstand-wert));
else
d=sqrt(dx*dx+dy*dy+(abstand)*(abstand));
if(inf.flag==1)
return(exp(-sqrphi_scal*d)/d);
else
return(wert*exp(-sqrphi_scal*d)/d);
#if ABSCHATTUNG
}
else{
flagfeld[j*inf.dim+i]=1; // abgeschattet
return 0.0;
}
#endif
}

// diese spezielle Version von getstrom erhöht außerdem die
// Variable anz_ring_ok, falls nicht abgeschattet wird, für
// checkdim
doub getstrom_special(int x, int y, int i, int j, int *anz){
doub d, wert=probfeld[j*inf.dim+i]; // wert= Funktionswert
int index; // am Aufpunkt
int dx=i-x;
int dy=j-y;
int dist=rint(sqrt(dx*dx+dy*dy));

index=floor((abstand-1.001*kor)*(doub)dist/(abstand-wert));
/* floor und "1.001" sollen nur einen
   Sicherheitsabstand gewährleisten */
while(index < dist &&
  abstand-(doub)index*(abstand-wert)/dist >
  LINIE_Z(x,y,dx,dy,dist,index))
index++;
if(dist==index){
++*anz; // keine Abschattung
flagfeld[j*inf.dim+i]=2; // 2: schon berechnet
//d=sqrt(dx*dx+dy*dy+(abstand)*(abstand));
d=sqrt(dx*dx+dy*dy+(abstand-wert)*(abstand-wert));
return (wert*exp(-sqrphi_scal*d)/d);
}
else{
flagfeld[j*inf.dim+i]=1; // der eigentliche
return 0.0; // Punkt selbst
}
}

// setzt das flagfeld auf Null
void clearflagfeld(){
register i;
int help=inf.dim*inf.dim;
for(i=0; i<help; i++)
flagfeld[i]=0;
}

// hier wird der Anteil der nicht abgeschatteten Punkte
// im äußersten Kreisring abgeschätzt
int checkdim(){
int diff; // halbe Sehnenlänge
int anz_ring; // Anzahl der Punkte in einem Ring
int anz_ring_ok; // Anzahl der Punkte ohne Abschattung
doub cur_ring; // Strom eines Kreisrings
doub cur=0; // gesamter Strom
int punkte=0; // gesamte Anzahl der Punkte
doub dr=(doub)inf.rmax/10.0; /* Delta-r: Schrittweite
   beim Vergrößern des
   Radius */
doub r=dr; /* Radius: mit dr wird begonnen,
   dann schrittweise erhöht */
doub r_alt; // der vorige Radius
int x=rint(inf.dim/2.0); // Spalte
int y=rint(inf.dim/2.0); // Zeile
// Es wird einfach der mittlere Punkt genommen !
int i, j;

clearflagfeld(); // Flags auf Null setzen

r_alt=0; /* das erste Gebiet ist kein Kreisring,
   sondern eine Kreisfläche */

while(r <= (doub)inf.rmax){
j=y-rint(r);
anz_ring=0;
anz_ring_ok=0; // es beginnt ein neuer Ring
cur_ring=0;
while(j<=y+rint(r)){
i=0; // es wird erstmal ganz links angefangen
// gehe bis an den linken Rand des Kreises:
while((i-x)*(i-x)+(j-y)*(j-y) > rint(r)*rint(r) &&
  ++i <= x);
if (0 > (diff=x-i)) // diff= Abstand von der y-Achse
exit(3); // man weiß ja nie
while(i<=x+diff){
if(r_alt && i<=x &&
   (i-x)*(i-x)+(j-y)*(j-y) <=
   rint(r_alt)*rint(r_alt))
i += 2*(x-i)+1; /* überspringe die
   innere Kreisfläche */
if (!flagfeld[j*inf.dim+i])
cur_ring+=getstrom_special(x, y, i, j,
   &anz_ring_ok);
anz_ring++; /* die Anzahl der Punkte
   im aktuellen Ring */
i++;
} // i-Schleife
j++;
} // j-Schleife
cur+=cur_ring;
punkte+=anz_ring;
r_alt=r;
r+=dr;
} // r-Schleife
inf.ant=(doub)anz_ring_ok / (doub)anz_ring; // nur zur Info
return 0; // Erfolg
}

// berechnet den Strom für den Punkt (x,y)
doub berechne(int x, int y){
register int i, j;
doub cur=0;
int diff;
clearflagfeld(); // Flags auf Null setzen

for(j=y-inf.rmax; j<=y+inf.rmax; j++){
i=0; // es wird erstmal ganz links angefangen
// gehe bis an den linken Rand des Kreises:
while((i-x)*(i-x)+(j-y)*(j-y) > inf.rmax*inf.rmax &&
  i++ <= x);
if (0 > (diff=x-i)) /* diff = Abstand von
   der y-Achse */
exit(4); // man weiß ja nie
while(i<=x+diff){
if (!flagfeld[j*inf.dim+i])
cur+=getstrom(x, y, i, j);
i++;
} // i-Schleife
} // j-Schleife
return cur;
}

// Hauptfuktion für die Rechnung
void rechne(char *kontrastfile){
register int i, j;

FILE *fp;
doub *feld;
doub max, min, hilf2;
char *tempfeld, *tempfeld2;
char str[DAT_NAME_MAX+20];
char str2[DAT_NAME_MAX+20];

abstand = (inf.c[1]/2+inf.d0)*inf.scal; // inf.d0 ist
// der Abstand
git = inf.scal*inf.c[0]; /* von Kern zu
   Spitze */
kor = inf.scal*inf.c[1];

hilf2 = log(1.0 - inf.cut/100)/sqrt(inf.phi);
inf.rmax= rint(inf.scal*sqrt(hilf2*hilf2 - 2*hilf2*inf.d0));
inf.dim= (int)rint(3*git + 2*inf.rmax +3);

feld =(doub *)malloc(SOD*inf.dim*inf.dim);
probfeld =(doub *)malloc(SOD*inf.dim*inf.dim);
flagfeld =(short *)calloc(inf.dim*inf.dim, SOS);
sqrphi_scal=sqrt(inf.phi)/inf.scal; // Hilfsvariable
sq3=sqrt(3); // Hilfsvariable

hilf=M_PI/git; // für das Makro surfunc

#if DEBUG
printf("inf.cut=%lf\ninf.rmax=%d\ninf.dim=%d\n",
   inf.cut, inf.rmax, inf.dim);
printf("abstand=%lf\ngit=%lf\nkor=%lf\n",
   abstand, git, kor);
#endif

surface(probfeld);

// checkdim schätzt inf.ant ab

checkdim();

for(j = 0; j < inf.dim; ++j)
for(i = 0; i < inf.dim; ++i)
feld[inf.dim*j+i] =
inf.rmax<=j && j<inf.dim-inf.rmax &&
inf.rmax<=i && i<inf.dim-inf.rmax ?
berechne(i, j) : probfeld[inf.dim*j + i];

/* ---------------------------------- */

/*   Es folgt das Erzeugen der Dateien xxx.tip,
  xxx.sur und xxx.sim */

/* noch ein letztes Mal updaten, bevor Benutzung! */
sprintf(inf.doc,
"%1.1lf#%1.2lf#%1.6lf#%d#%1.1lf#%1.1lf#%1.1lf#%1.1lf",
inf.phi,inf.ant,inf.cut,inf.rmax,inf.scal,inf.c[0],inf.c[1],
inf.d0);

// xxx.sur
j=inf.dim*inf.dim;
tempfeld=tempfeld2=(char *)malloc(SOC*j);
max=maximum(probfeld, j);
min=minimum(probfeld, j);
for(i=0; i<j; i++)
tempfeld[i]=(int)((probfeld[i]-min)/(max-min)*255);
sprintf(str, "%s.sur", inf.doc);
#if 0
if(!(fp=fopen(str, "w")))
exit(1);
schreibedim(inf.dim, fp);
fwrite(tempfeld2, SOC, j, fp);
#endif

// xxx.sim
tempfeld=(char *)malloc(SOC*j);
max=feld[inf.rmax*inf.dim+inf.rmax];
min=feld[inf.rmax*inf.dim+inf.rmax];
for(j=inf.rmax; j<inf.dim-inf.rmax; j++)
for(i=inf.rmax; i<inf.dim-inf.rmax; i++){
if(feld[j*inf.dim+i] > max)
max=feld[j*inf.dim+i];
if(feld[j*inf.dim+i] < min)
min=feld[j*inf.dim+i];
}
#if DEBUG
printf("\n---> max=%le und min=%le <---\n", max, min);
printf("\n---> inf.rmax=%d und inf.rmax=%d <---\n",
   inf.rmax, inf.rmax);
#endif
for(j=0; j<inf.dim; j++)
for(i=0; i<inf.dim; i++)
tempfeld[j*inf.dim+i]= inf.rmax<=j &&
j<inf.dim-inf.rmax && inf.rmax<=i &&
i<inf.dim-inf.rmax ?
(int)((feld[j*inf.dim+i]-min)/
  (max-min)*255) :
tempfeld2[j*inf.dim+i];
sprintf(str, "%s.sim", inf.doc);
fp=fopen(str, "w");
schreibedim(inf.dim, fp); // provisorisch
fwrite(tempfeld, SOC, j*j, fp);

// Erst wird der Name für die Kontrastdatei erzeugt:
// (falls keine Datei benannt wurde wird "kontrast" benutzt)
#if ABSCHATTUNG
sprintf(str, "%s", *kontrastfile ? kontrastfile : "kont");
#else
sprintf(str, "%s", *kontrastfile ? kontrastfile : "kontoa");
#endif
if(!(*kontrastfile)){
sprintf(str2, ".%2.0lf", inf.cut);
strcat(str, str2);
}

freopen(str, "a", fp);
fprintf(fp, "%s: %lf %lf\n", inf.doc, inf.d0, 1.0- min/max);
fclose(fp);
#if HP
sprintf(str, "gzip %s.sim &", inf.doc);
system(str);
#endif

// Speicher freigeben:
free(feld);
free(flagfeld);
free(probfeld);
free(tempfeld);
free(tempfeld2);
}

/* simserie.c

Author: Peter Münster, 1996

Enthält Funktionen für Rechenserien

*/

#include "sim.h"

extern info_struct inf;
extern char progname[];

// paßt Parameter an und startet rechne

void serienfunc(char *str){
// inf.rmax=inf.d0 * tan(M_PI*inf.wmax /180);
// inf.dim = 3*(int)maximum(inf.c, 4) +
// 3+2*(int)(inf.scal*inf.rmax);
/* Das gescannte Feld beträgt also 3 Gitterkonstanten */
/* Der Summand dims ist nötig um die
Einflüsse von außerhalb */
/* des gescannten Feldes mit zu berücksichtigen */
rechne(str);
}

// vorläufig wird nur der Abstand variiert!

void serie(proc_struct **proc_ptr, int *count){
char str[DAT_NAME_MAX];
pid_t fork_ret;

my_doub_scan("\nKleinster Abstand?\n[%lf]: ", &inf.d0);
my_doub_scan("\nGrößter Abstand?\n[%lf]: ", &inf.d0max);
my_doub_scan("\nSchrittweite?\n[%lf]: ", &inf.dd0);

pf1(
"\nNamen der Kontrastdatei eingeben!\n(default=%s): ",
"kontrast");
if(*((fgets(str, DAT_NAME_MAX, stdin))+1)) // wegen Return
str[strlen(str)-1]=’\0’;
else
*str=’\0’;

fork_ret=fork();
switch(fork_ret){
case -1:
FPF("fork-Aufruf in serie %s", "ohne Erfolg!\n");
break;
case 0:
for(; inf.d0<inf.d0max; inf.d0+=inf.dd0)
serienfunc(str);
exit(0);
default:
(*count)++; /* Prozeßinfos abspeichern: */
*proc_ptr=(proc_struct *)realloc(*proc_ptr,
sizeof(proc_struct)* *count);
(*proc_ptr)[*count-1].pid=fork_ret;
strncpy((*proc_ptr)[*count-1].name,
str, DAT_NAME_MAX);
(*proc_ptr)[*count-1].name[DAT_NAME_MAX-1]=’\0’;
}
}

/* simhilf.c

Author: Peter Münster, 1996

Hilfsfunktionen für sim

*/

#include "sim.h"

extern char *progname;

#if 0
// getch() mit echo

int getche()
{
int c;
noecho();
cbreak();
c=getch();
nocbreak();
echo();
if (c == 127)
mywrite("\b");
else
pf1("%c", c);
return(c);
}

// getch() ohne echo

int getchne()
{
int c;
noecho();
cbreak();
c=getch();
nocbreak();
echo();
return(c);
}
#endif

// von Max

void schreibedim(int dim, FILE *ausdp){
/* schreibedim schreibt die Dimension zweimal */
/* d.h. es sind quadratische Bilder */
int i, low, high;

low = dim%256;
high = (dim - low)/256;
for(i = 0; i < 2; ++i){
putc(high,ausdp);
putc(low,ausdp);
}
}

// maximum des valeurs de p

doub maximum(doub *p, int j){
register unsgnd i;
doub m;

m = *p;
for (i=0; i<j; ++i)
if (m < *++p) m = *p;
return m;
}

// minimum des valeurs de p

doub minimum(doub *p, int j){
register unsgnd i;
doub m;

m = *p;
for (i=0; i<j; ++i)
if (m > *++p) m = *p;
return m;
}

Mit folgender kleiner Beispiel-Makedatei kann das Programm bequem übersetzt werden:

B.1. Simulationsprogramme zur Kontrastinversion