1
2 /****************************************************************************
3 **
4 ** Copyright (C) 2011 Christian B. Huebschle & George M. Sheldrick
5 ** All rights reserved.
6 ** Contact: chuebsch@moliso.de
7 **
8 ** This file is part of the ShelXle
9 **
10 ** This file may be used under the terms of the GNU Lesser
11 ** General Public License version 2.1 as published by the Free Software
12 ** Foundation and appearing in the file COPYING included in the
13 ** packaging of this file. Please review the following information to
14 ** ensure the GNU Lesser General Public License version 2.1 requirements
15 ** will be met: http://www.gnu.org/licenses/old-licenses/lgpl-2.1.html.
16 **
17 **
18 ****************************************************************************/
19 #include "fourxle.h"
20 #include <QtCore>
21 #include <QtGui>
22 #include "deprecation.h"
23 #include "itsme.h"
24 //üäöß
FourXle(Molecule * mole_,ChGL * chgl_,QToolBar * toolView,double resol,double wght)25 FourXle::FourXle(Molecule *mole_, ChGL *chgl_,QToolBar *toolView, double resol, double wght){
26 // printf("FourXle\n");
27 ownfcf6_=nullptr;
28 mole =mole_;
29 f000=0.0;
30 overridef000=-1.0;
31 n1=n2=n3=n4=n5=0;
32 voxelstr.clear();
33 chgl=chgl_;
34 MYFCF=false;
35 oldatomsize=-1;
36 chgl->habzutun=false;
37 map_radius=5.0;
38 maptrunc = 2;
39 mapcontrol=nullptr;
40 resDir=QDir::current();
41 // for (int i=0; i<37; i++) chgl->foubas[i]=0;
42 /*
43 QVector< QVector< float > > chgl->fchgl->fVertexes;
44 QVector< QVector< float > > chgl->fNormals;
45 */
46 chgl->fVertexes.resize(28);
47 chgl->fNormals.resize(28);
48 HKLMX=200;
49 datfo_fc=datfo=nullptr;
50 nodex=nodey=nodez=nullptr;noGoMap=nullptr;
51 nodeX=nodeY=nodeZ=nullptr;
52 urs=V3(0,0,0);
53 nr=0;
54 nc=0;
55 sigma[0]=sigma[1]=iso[0]=iso[1]=0;
56 chgl->lintrans=0.8;
57 chgl->linwidth=0.5;
58 rr=resol;
59 rw=wght;
60 chgl->foact=toolView->addAction(QIcon(":/xle-icons/fomaps.svg"),"toggle Fo map",chgl,SLOT(updateGL()));
61 chgl->foact->setWhatsThis("<h1>Toggle Fo map <img src=\":/xle-icons/fomaps.svg\"></h1> Toggle F-observed map on/off.");
62 connect(chgl->foact,SIGNAL(destroyed(QObject*)),this,SLOT(foactDestroyed()));
63 chgl->foact->setCheckable(true);
64 chgl->fofcact=toolView->addAction(QIcon(":/xle-icons/diffmap.svg"),"toggle Fo-Fc map",chgl,SLOT(updateGL()));
65 chgl->fofcact->setWhatsThis("<h1>Toggle Fo-Fc map<img src=\":/xle-icons/diffmap.svg\"></h1> Toggle (F-observed)-(F-calculated) map on/off. Fo-Fc map is also known as difference Fourier map.");
66 connect(chgl->fofcact,SIGNAL(destroyed(QObject*)),this,SLOT(fofcactDestroyed()));
67 chgl->fofcact->setCheckable(true);
68 chgl->fofcact->setVisible(false);
69 chgl->foact->setVisible(false);
70 doMaps = new QCheckBox("Calculate Maps");
71 doMaps->setChecked(true);
72 doMaps->setWhatsThis("<h1>Calculate maps</h1> Enables Fo and Fo-Fc map calculation. This can only be performed if LIST is set to 6.");
73 keepIso = new QCheckBox("Keep map iso values");
74 keepIso->setChecked(false);
75 keepIso->setVisible(false);
76 connect(keepIso,SIGNAL(toggled(bool)), this,SLOT(keepIsoValues(bool )));
77 connect(chgl,SIGNAL(destroyed(QObject*)),this,SLOT(chglDestroyed()));
78 char *ENUS=setlocale(LC_ALL,"C");
79 if (!ENUS)ENUS=setlocale(LC_ALL,"en_US.UTF-8");
80 {
81 buttonBoxMC = new QDialogButtonBox( QDialogButtonBox::Apply | QDialogButtonBox::Cancel);
82 applyMC = buttonBoxMC->button(QDialogButtonBox::Apply);
83 applyDF = buttonBoxMC->addButton("Apply defaults",QDialogButtonBox::ActionRole);
84 wl= new QLabel("Factor to downweight weak data");
85 pl= new QLabel("Map precision: ");
86 dl= new QLabel("Fo-Fc map: ");
87 ol= new QLabel("F-observed map: ");
88 sl= new QLabel("Map truncation type: ");
89 lw= new QLabel("Line width: ");
90 ltz= new QLabel("Line transparency: ");
91 md = new QDialog();
92 md->setWindowTitle("Map Control");
93 connect(applyMC , SIGNAL(clicked()), this , SLOT(controlMap()));
94 connect(applyMC , SIGNAL(clicked()), this , SLOT(openMapControl()));
95 connect(applyDF , SIGNAL(clicked()), this , SLOT(mapDefault()));
96 connect(buttonBoxMC, SIGNAL(rejected()), md, SLOT(hide()));
97
98 mapprec = new QDoubleSpinBox(md);
99 mapprec->setMinimum(1.1);
100 mapprec->setMaximum(8.0);
101 mapprec->setValue(rr);
102 mapprec->setSingleStep(0.1);
103
104 lineTrans = new QDoubleSpinBox(md);
105 lineTrans->setMinimum(0.1);
106 lineTrans->setMaximum(1.0);
107 lineTrans->setValue(chgl->lintrans);
108 lineTrans->setSingleStep(0.1);
109
110 lineWidth = new QDoubleSpinBox(md);
111 lineWidth->setMinimum(0.15);
112 lineWidth->setMaximum(4.0);
113 lineWidth->setValue(chgl->linwidth);
114 lineWidth->setSingleStep(0.1);
115
116 weak = new QDoubleSpinBox(md);
117 weak->setMinimum(0.1);
118 weak->setMaximum(8.0);
119 weak->setValue(rw);
120 weak->setSingleStep(0.1);
121
122 difmaps = new QDoubleSpinBox(md);
123 difmaps->setMinimum(0.025);
124 difmaps->setMaximum(9999);
125 difmaps->setValue(fabs(iso[1]));
126 difmaps->setSingleStep(fabs(iso[1])/10.0);
127 difmaps->setSuffix("e/A^3");
128
129 fomaps = new QDoubleSpinBox();
130 fomaps->setMinimum(0.025);
131 fomaps->setMaximum(9999);
132 fomaps->setValue(iso[0]);
133 fomaps->setSingleStep(iso[0]/10.0);
134 fomaps->setSuffix("e/A^3");
135
136 cutrange = new QDoubleSpinBox(md);
137 cutrange->setMinimum(0.5);
138 cutrange->setMaximum(3.0);
139 cutrange->setValue(1.41);
140 cutrange->setSingleStep(0.1);
141 cutrange->setSuffix("A");
142 cutrange->setPrefix("cut range ");
143
144 maprad = new QDoubleSpinBox();
145 maprad->setValue(map_radius);
146 maprad->setMinimum(0);
147 maprad->setMaximum(25.0);
148 maprad->setSuffix("A");
149 maprad->setPrefix("Map radius: ");
150 maprad->setSingleStep(0.3);
151 mapSchnitt= new QComboBox();
152 mapSchnitt->addItem("one complete unit cell",0);
153 mapSchnitt->addItem("sphere around rotation center",1);
154 mapSchnitt->addItem("cut range around visible atoms or peaks",2);
155 mapSchnitt->setCurrentIndex(maptrunc);
156 connect(mapSchnitt,SIGNAL(currentIndexChanged(int)),this,SLOT(schnittart(int)));
157 chgl->lighting = new QCheckBox("Lighting");
158 chgl->lighting->setChecked(true);
159 chgl->niceTrans= new QCheckBox("Undistorted transparency");
160 chgl->niceTrans->setChecked(false);
161 chgl->fillMap= new QCheckBox("Filled map surface");
162 chgl->fillMap->setChecked(false);
163 mdl = new QGridLayout();
164 mdl->addWidget(mapprec,0,1);
165 mdl->addWidget(weak,1,1);
166 mdl->addWidget(difmaps,2,1);
167 mdl->addWidget(fomaps,3,1);
168 mdl->addWidget(maprad,6,2);
169 mdl->addWidget(mapSchnitt,6,1);
170 mdl->addWidget(cutrange,6,2);
171 mdl->addWidget(lineTrans,7,1);
172 mdl->addWidget(lineWidth,8,1);
173 mdl->addWidget(chgl->lighting,9,1);
174 mdl->addWidget(chgl->niceTrans,9,0);
175 mdl->addWidget(chgl->fillMap,9,2);
176 connect(lineTrans,SIGNAL(valueChanged(double)),this,SLOT(updateLines()));
177 connect(lineWidth,SIGNAL(valueChanged(double)),this,SLOT(updateLines()));
178
179 connect(chgl->fillMap,SIGNAL(stateChanged(int)),chgl,SLOT(updateGL()));
180 connect(chgl->niceTrans,SIGNAL(stateChanged(int)),chgl,SLOT(updateGL()));
181 connect(chgl->lighting,SIGNAL(stateChanged(int)),chgl,SLOT(updateGL()));
182 mdl->addWidget(pl,0,0);
183 mdl->addWidget(wl,1,0);
184 mdl->addWidget(dl,2,0);
185 mdl->addWidget(ol,3,0);
186 //mdl->addWidget(rl,5,0);
187 mdl->addWidget(sl,6,0);
188 mdl->addWidget(ltz,7,0);
189 mdl->addWidget(lw,8,0);
190
191 QTabWidget *tabWidget = new QTabWidget(md);
192 QWidget *gwmd,*admd;
193 gwmd = new QWidget();
194
195 gwmd->setLayout(mdl);
196 admd = new QWidget();
197 QVBoxLayout *mdmain= new QVBoxLayout;
198 tabWidget->addTab(gwmd,"General");
199 tabWidget->addTab(admd,"Advanced");
200
201 mole->einstellung->beginGroup("Colors");
202 chgl->fopc=mole->einstellung->value("FourFoPlusColor", QColor("#0000FF")).value<QColor>();
203 chgl->fomc=mole->einstellung->value("FourFoMinusColor", QColor("#FFA600")).value<QColor>();
204 chgl->dipc=mole->einstellung->value("FourDiffPlusColor", QColor("#00B200")).value<QColor>();
205 chgl->dimc=mole->einstellung->value("FourDiffMinusColor",QColor("#CC0000")).value<QColor>();
206 mole->einstellung->endGroup();
207
208 foPlusButton =new QPushButton("Fo positive color");
209 foPlusButton->setStyleSheet(QString(
210 "QPushButton {"
211 "font-size:16px;"
212 "font-weight:bold;"
213 "border: 1px solid #000000;"
214 "border-radius: 9px;"
215 "color: %4;"
216 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
217 "}"
218 "QPushButton:hover {"
219 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
220 "border: 0px"
221 "}"
222 "QPushButton:flat {"
223 " border: none; /* no border for a flat push button */"
224 "}"
225 )
226 .arg( chgl->fopc.name())
227 .arg( chgl->fopc.darker(200).name())
228 .arg( chgl->fopc.lighter(100).name())
229 .arg("#ffffff"));
230 connect(foPlusButton,SIGNAL(pressed()),this,SLOT(colorDLGFOP()));
231 defColorButton =new QPushButton("Default color scheme");
232 connect(defColorButton,SIGNAL(pressed()),this,SLOT(defaultColors()));
233 foMinusButton =new QPushButton("Fo negative color");
234 foMinusButton->setStyleSheet(QString(
235 "QPushButton {"
236 "font-size:16px;"
237 "font-weight:bold;"
238 "border: 1px solid #000000;"
239 "border-radius: 9px;"
240 "color: %4;"
241 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
242 "}"
243 "QPushButton:hover {"
244 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
245 "border: 0px"
246 "}"
247 "QPushButton:flat {"
248 " border: none; /* no border for a flat push button */"
249 "}"
250 )
251 .arg( chgl->fomc.name())
252 .arg( chgl->fomc.darker(200).name())
253 .arg( chgl->fomc.lighter(100).name())
254 .arg("#ffffff"));
255 connect(foMinusButton,SIGNAL(pressed()),this,SLOT(colorDLGFOM()));
256
257 diffMinusButton =new QPushButton("Fo-Fc negative color");
258 diffMinusButton->setStyleSheet(QString(
259 "QPushButton {"
260 "font-size:16px;"
261 "font-weight:bold;"
262 "border: 1px solid #000000;"
263 "border-radius: 9px;"
264 "color: %4;"
265 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
266 "}"
267 "QPushButton:hover {"
268 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
269 "border: 0px"
270 "}"
271 "QPushButton:flat {"
272 " border: none; /* no border for a flat push button */"
273 "}"
274 )
275 .arg( chgl->dimc.name())
276 .arg( chgl->dimc.darker(200).name())
277 .arg( chgl->dimc.lighter(100).name())
278 .arg("#ffffff"));
279 connect(diffMinusButton,SIGNAL(pressed()),this,SLOT(colorDLGDIM()));
280
281 diffPlusButton =new QPushButton("Fo-Fc positive color");
282 diffPlusButton->setStyleSheet(QString(
283 "QPushButton {"
284 "font-size:16px;"
285 "font-weight:bold;"
286 "border: 1px solid #000000;"
287 "border-radius: 9px;"
288 "color: %4;"
289 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
290 "}"
291 "QPushButton:hover {"
292 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
293 "border: 0px"
294 "}"
295 "QPushButton:flat {"
296 " border: none; /* no border for a flat push button */"
297 "}"
298 )
299 .arg( chgl->dipc.name())
300 .arg( chgl->dipc.darker(200).name())
301 .arg( chgl->dipc.lighter(100).name())
302 .arg("#ffffff"));
303 connect(diffPlusButton,SIGNAL(pressed()),this,SLOT(colorDLGDIP()));
304
305 amdl = new QVBoxLayout();
306 amdl->addWidget(foPlusButton);
307 amdl->addWidget(foMinusButton);
308 amdl->addWidget(diffPlusButton);
309 amdl->addWidget(diffMinusButton);
310 amdl->addWidget(defColorButton);
311 jnkButton = new QPushButton("fractal dimension analysis");
312 amdl->addWidget(jnkButton);
313 admd->setLayout(amdl);
314 mdmain->addWidget(tabWidget);
315 mdmain->addWidget(buttonBoxMC);
316
317 md->setLayout(mdmain);
318 md->setModal (false);
319 ownfcf6=new QCheckBox("fcf6 by hkl");
320 ownfcf6->setChecked(MYFCF);
321 ownfcf6_=new QCheckBox("Compute fcf6 file from structure factors and observed data from hkl file");
322 ownfcf6_->setChecked(MYFCF);
323 amdl->addWidget(ownfcf6_);
324 connect(ownfcf6_,SIGNAL(toggled(bool)),this , SLOT(doMapsNow(bool)));
325 connect(doMaps,SIGNAL(toggled(bool)),this , SLOT(doMapsNow(bool)));
326 mole->einstellung->beginGroup("Window");//mole->einstellung->beginGroup("");
327 if (mole->einstellung->contains("Fo-Map")) chgl->foact->setChecked(mole->einstellung->value("Fo-Map").toBool());
328 if (mole->einstellung->contains("Fc-Fo-Map")) chgl->fofcact->setChecked(mole->einstellung->value("Fc-Fo-Map").toBool());
329 if (mole->einstellung->contains("Map-truncation-type")) maptrunc=mole->einstellung->value("Map-truncation-type").toInt();
330 if (mole->einstellung->contains("Map-Line-Transparency")) chgl->lintrans = mole->einstellung->value("Map-Line-Transparency").toDouble();
331 if (mole->einstellung->contains("Map-Line-Width")) chgl->linwidth = mole->einstellung->value("Map-Line-Width").toDouble();
332 if (mole->einstellung->contains("Map-Resolution")) rr = mole->einstellung->value("Map-Resolution").toDouble();
333 if (mole->einstellung->contains("Map-Weight")) rw = mole->einstellung->value("Map-Weight").toDouble();
334 if (mole->einstellung->contains("Map-Lighting")) chgl->lighting->setChecked(mole->einstellung->value("Map-Lighting").toBool());
335 if (mole->einstellung->contains("Nice-Transparency")) chgl->niceTrans->setChecked(mole->einstellung->value("Nice-Transparency").toBool());
336 if (mole->einstellung->contains("Solid-Map")) chgl->fillMap->setChecked(mole->einstellung->value("Solid-Map").toBool());
337 if (mole->einstellung->contains("Calculate-Maps")) doMaps->setChecked(mole->einstellung->value("Calculate-Maps").toBool());
338 if (mole->einstellung->contains("KeepMapIsoValues")){
339 keepIso->setChecked(mole->einstellung->value("KeepMapIsoValues").toBool());
340 // keepiso->setChecked(mole->einstellung->value("KeepMapIsoValues").toBool());
341 }
342 mole->einstellung->endGroup();
343 schnittart(maptrunc);
344 // printf("@!#$ %p\n",mole->einstellung);
345 md->hide();
346 }
347 connect(ownfcf6,SIGNAL(toggled(bool)),this , SLOT(doMapsNow(bool)));
348 isBelo=false;
349 //printf("constr. belo?%d\n",isBelo);
350 //printf("%s\n",ENUS);
351 }
352
~FourXle()353 FourXle::~FourXle(){
354 //printf("Maps will close...\n");
355 if (chgl!=nullptr) disconnect(chgl,SIGNAL(diffscroll(int ,int )),0,0);
356 if (chgl!=nullptr) disconnect(chgl,SIGNAL(neuemitte(V3)),0,0);
357 if (chgl!=nullptr) disconnect(chgl,SIGNAL(inimibas()),0,0);
358 //printf("%d\n",__LINE__);
359 killmaps();
360 //printf("%d\n",__LINE__);
361 delete doMaps;
362 //printf("Maps closed.\n");
363 }
364
mapDefault()365 void FourXle::mapDefault(){
366 map_radius= 5.0;
367 maptrunc = 2;
368 rr= 2.0;
369 rw= 1.0;
370 chgl->lintrans= 0.8;
371 chgl->linwidth= 0.5;
372 iso[1]= sigma[0] * 2.7f;
373 iso[0]= sigma[1] * 1.2f;
374 chgl->niceTrans->setChecked(false);
375 chgl->lighting->setChecked(true);
376 chgl->fillMap->setChecked(false);
377 weak->setValue(rw);
378 mapprec->setValue(rr);
379 difmaps->setValue(fabs(iso[1]));
380 difmaps->setSingleStep(fabs(iso[1])/10.0);
381 lineTrans->setValue(chgl->lintrans);
382 lineWidth->setValue(chgl->linwidth);
383 fomaps->setValue(iso[0]);
384 fomaps->setSingleStep(iso[0]/10.0);
385 maprad->setValue(map_radius);
386 mapSchnitt->setCurrentIndex(maptrunc);
387 controlMap();
388 }
389
schnittart(int trunc)390 void FourXle::schnittart(int trunc){
391 if (trunc==1) maprad->show();
392 else maprad->hide();
393 if (trunc==2) cutrange->show();
394 else cutrange->hide();
395 }
396
keepIsoValues(bool b)397 void FourXle::keepIsoValues(bool b){
398 keepiso->setChecked(b);
399 keepIso->setChecked(b);
400 keepIso->setVisible(b);
401 if (n5>0){
402 sigma[0] = sigmaFoFc();
403 sigma[1] = sigmaFo();
404 iso[1]= sigma[0] * 2.7f;
405 iso[0]= sigma[1] * 1.2f;
406 fomaps->setValue(iso[0]);
407 difmaps->setValue(fabs(iso[1]));
408 controlMap();
409 }
410 }
411
updateLines()412 void FourXle::updateLines(){
413 chgl->lintrans=lineTrans->value();
414 chgl->linwidth=lineWidth->value();
415 chgl->updateGL();
416 }
417
killFourXleFirst()418 void FourXle::killFourXleFirst(){
419 //fprintf(stderr,"chgl dies...argh...\n");
420 killmaps();
421 n5=0;
422 if (mapcontrol!=nullptr) mapcontrol->setVisible(false);
423 }
424
doMapsNow(bool b)425 void FourXle::doMapsNow(bool b){
426 //isBelo=false;
427 bool bb=MYFCF;
428 bool B = (MYFCF!=ownfcf6->isChecked())||(MYFCF!=ownfcf6_->isChecked());
429 MYFCF=(B)?!MYFCF:MYFCF;
430 if (bb!=MYFCF) {
431 ownfcf6->disconnect();
432 ownfcf6_->disconnect();
433 ownfcf6->setChecked(MYFCF);
434 ownfcf6_->setChecked(MYFCF);
435 connect(ownfcf6,SIGNAL(toggled(bool)),this , SLOT(doMapsNow(bool)));
436 connect(ownfcf6_,SIGNAL(toggled(bool)),this , SLOT(doMapsNow(bool)));
437 emit reloadRes();
438 return;
439 }
440 if (!b) {
441 killmaps();
442 n5=0;
443 if (mapcontrol!=nullptr) mapcontrol->setVisible(false);
444 chgl->updateGL();
445 return;
446 } else {
447 bool hamrsho=((iso[0]==iso[1])&&((iso[0]==0.0f)||(iso[0]==0.03f)));
448 if (!loadFouAndPerform(fouName.toLocal8Bit(),hamrsho)){
449 emit bigmessage(QString("<font color=red>Error while loading %1!</font><br>").arg(fouName));
450 if (mapcontrol!=nullptr) mapcontrol->setVisible(false);
451 } else {
452 if (mapcontrol!=nullptr) mapcontrol->setVisible(true);
453 if (!hamrsho) {
454 sigma[0] = sigmaFoFc();
455 sigma[1] = sigmaFo();
456 iso[1]= sigma[0] * 2.7f;
457 iso[0]= sigma[1] * 1.2f;
458 fomaps->setValue(iso[0]);
459 difmaps->setValue(fabs(iso[1]));
460 inimap();
461 }
462 }
463 chgl->updateGL();
464 }
465
466 }
467
controlMap()468 void FourXle::controlMap(){
469 maptrunc=mapSchnitt->currentIndex ();
470 if ((rw!=weak->value())||(rr!=mapprec->value())){
471 rr = mapprec->value();
472 rw = weak->value();
473 if (!loadFouAndPerform(fouName.toLocal8Bit(),false)){
474 emit bigmessage(QString("<font color=red>Error while loading %1!</font><br>").arg(fouName));
475 }
476 }
477 chgl->lintrans = lineTrans->value();
478 chgl->linwidth = lineWidth->value();
479 map_radius = maprad->value();
480 iso[0] = fomaps->value();
481 iso[1] = -difmaps->value();
482 inimap();
483 chgl->updateGL();
484 }
485
killmaps()486 void FourXle::killmaps(){
487 //printf("KM belo?%d chgl->%p dataFo-Fc->%p dataFo->%p\n",isBelo,chgl,datfo_fc,datfo);
488 isBelo=false;
489 // printf("These nodes will be killed %p %p %p\n",nodex,nodey,nodez);
490 /*! deletes all fourier maps and frees the memory
491 */
492 // printf("I kill maps\n");
493
494 if (chgl!=nullptr) disconnect(chgl,SIGNAL(diffscroll(int ,int )),0,0);
495 if (chgl!=nullptr) disconnect(chgl,SIGNAL(neuemitte(V3)),0,0);
496 if (chgl!=nullptr) disconnect(chgl,SIGNAL(inimibas()),0,0);
497 //printf("emit halt\n");
498 //emit haltMW();
499 //printf("halt emited\n");
500 if ((chgl!=nullptr)&&(chgl->fofcact!=nullptr)) chgl->fofcact->setVisible(false);
501 if ((chgl!=nullptr)&&(chgl->foact!=nullptr)) chgl->foact->setVisible(false);
502 if (datfo!=nullptr) free(datfo);
503 if (datfo_fc!=nullptr) free(datfo_fc);
504 if (nodex!=nullptr) free(nodex);
505 nodex=nullptr;
506 if (nodey!=nullptr) free(nodey);
507 nodey=nullptr;
508 if (nodez!=nullptr) free(nodez);
509 nodez=nullptr;
510 if (noGoMap!=nullptr) free(noGoMap);
511 noGoMap=nullptr;
512 if (nodeX!=nullptr) free(nodeX);
513 if (nodeY!=nullptr) free(nodeY);
514 if (nodeZ!=nullptr) free(nodeZ);
515 nodeX=nodeY=nodeZ=nullptr;
516 datfo=datfo_fc=nullptr;
517 if (chgl!=nullptr){
518 //printf("I kill maps %d! %d!\n",chgl->fVertexes.size(),chgl->fNormals.size());
519 int sumsize=0;
520 for (int i=0; i<chgl->fVertexes.size(); i++) {sumsize+=chgl->fVertexes[i].size();chgl->fVertexes[i].clear();}
521 for (int i=0; i<chgl->fNormals.size(); i++) {sumsize+=chgl->fNormals[i].size();chgl->fNormals[i].clear();}
522 //printf("Freed %d normals and vertices.\n",sumsize);
523 }//else printf("chgl is already dead!\n");
524 }
525
colorDLGFOP()526 void FourXle::colorDLGFOP(){
527 QColor color;
528 color=QColorDialog::getColor(chgl->fopc, md);
529 if (color.isValid()) {
530 chgl->fopc=color;
531 foPlusButton->setStyleSheet(QString(
532 "QPushButton {"
533 "font-size:16px;"
534 "font-weight:bold;"
535 "border: 1px solid #000000;"
536 "border-radius: 9px;"
537 "color: %4;"
538 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
539 "}"
540 "QPushButton:hover {"
541 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %2, stop: 0.5 %1, stop: 1 %2);"
542 "}"
543 "QPushButton:flat {"
544 " border: none; /* no border for a flat push button */"
545 "}"
546 )
547 .arg(color.name())
548 .arg(color.darker(200).name())
549 .arg(color.lighter(100).name())
550 .arg("#ffffff"));
551 }
552 mole->einstellung->beginGroup("Colors");
553 mole->einstellung->setValue("FourFoPlusColor", chgl->fopc);
554 mole->einstellung->setValue("FourFoMinusColor", chgl->fomc);
555 mole->einstellung->setValue("FourDiffPlusColor", chgl->dipc);
556 mole->einstellung->setValue("FourDiffMinusColor",chgl->dimc);
557 mole->einstellung->endGroup();
558 md->update();
559
560 }
561
colorDLGFOM()562 void FourXle::colorDLGFOM(){
563 QColor color;
564 color=QColorDialog::getColor(chgl->fopc, md);
565 if (color.isValid()) {
566 chgl->fomc=color;
567 foMinusButton->setStyleSheet(QString(
568 "QPushButton {"
569 "font-size:16px;"
570 "font-weight:bold;"
571 "border: 1px solid #000000;"
572 "border-radius: 9px;"
573 "color: %4;"
574 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
575 "}"
576 "QPushButton:hover {"
577 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %2, stop: 0.5 %1, stop: 1 %2);"
578 "}"
579 "QPushButton:flat {"
580 " border: none; /* no border for a flat push button */"
581 "}"
582 )
583 .arg(color.name())
584 .arg(color.darker(200).name())
585 .arg(color.lighter(100).name())
586 .arg("#ffffff"));
587 }
588 mole->einstellung->beginGroup("Colors");
589 mole->einstellung->setValue("FourFoPlusColor", chgl->fopc);
590 mole->einstellung->setValue("FourFoMinusColor", chgl->fomc);
591 mole->einstellung->setValue("FourDiffPlusColor", chgl->dipc);
592 mole->einstellung->setValue("FourDiffMinusColor",chgl->dimc);
593 mole->einstellung->endGroup();
594 md->update();
595
596 }
597
colorDLGDIP()598 void FourXle::colorDLGDIP(){
599 QColor color;
600 color=QColorDialog::getColor(chgl->fopc, md);
601 if (color.isValid()) {
602 chgl->dipc=color;
603
604 diffPlusButton->setStyleSheet(QString(
605 "QPushButton {"
606 "font-size:16px;"
607 "font-weight:bold;"
608 "border: 1px solid #000000;"
609 "border-radius: 9px;"
610 "color: %4;"
611 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
612 "}"
613 "QPushButton:hover {"
614 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %2, stop: 0.5 %1, stop: 1 %2);"
615 "}"
616 "QPushButton:flat {"
617 " border: none; /* no border for a flat push button */"
618 "}"
619 )
620 .arg(color.name())
621 .arg(color.darker(200).name())
622 .arg(color.lighter(100).name())
623 .arg("#ffffff"));
624 }
625 mole->einstellung->beginGroup("Colors");
626 mole->einstellung->setValue("FourFoPlusColor", chgl->fopc);
627 mole->einstellung->setValue("FourFoMinusColor", chgl->fomc);
628 mole->einstellung->setValue("FourDiffPlusColor", chgl->dipc);
629 mole->einstellung->setValue("FourDiffMinusColor",chgl->dimc);
630 mole->einstellung->endGroup();
631 md->update();
632 }
633
colorDLGDIM()634 void FourXle::colorDLGDIM(){
635 QColor color;
636 color=QColorDialog::getColor(chgl->fopc, md);
637 if (color.isValid()) {
638 chgl->dimc=color;
639
640 diffMinusButton->setStyleSheet(QString(
641 "QPushButton {"
642 "font-size:16px;"
643 "font-weight:bold;"
644 "border: 1px solid #000000;"
645 "border-radius: 9px;"
646 "color: %4;"
647 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
648 "}"
649 "QPushButton:hover {"
650 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %2, stop: 0.5 %1, stop: 1 %2);"
651 "}"
652 "QPushButton:flat {"
653 " border: none; /* no border for a flat push button */"
654 "}"
655 )
656 .arg(color.name())
657 .arg(color.darker(200).name())
658 .arg(color.lighter(100).name())
659 .arg("#ffffff"));
660 }
661 mole->einstellung->beginGroup("Colors");
662 mole->einstellung->setValue("FourFoPlusColor", chgl->fopc);
663 mole->einstellung->setValue("FourFoMinusColor", chgl->fomc);
664 mole->einstellung->setValue("FourDiffPlusColor", chgl->dipc);
665 mole->einstellung->setValue("FourDiffMinusColor",chgl->dimc);
666 mole->einstellung->endGroup();
667 md->update();
668
669 }
670
openMapControl()671 void FourXle::openMapControl(){
672 weak->setValue(rw);
673 mapprec->setValue(rr);
674 difmaps->setValue(fabs(iso[1]));
675 difmaps->setSingleStep(fabs(iso[1])/10.0);
676 lineTrans->setValue(chgl->lintrans);
677 lineWidth->setValue(chgl->linwidth);
678 fomaps->setValue(iso[0]);
679 fomaps->setSingleStep(iso[0]/10.0);
680 maprad->setValue(map_radius);
681 mapSchnitt->setCurrentIndex(maptrunc);
682 md->show();
683 }
684
defaultColors()685 void FourXle::defaultColors(){
686 chgl->fopc = QColor("#0000FF");
687 chgl->fomc = QColor("#FFA600");
688 chgl->dipc = QColor("#00B200");
689 chgl->dimc = QColor("#CC0000");
690 foPlusButton->setStyleSheet(QString(
691 "QPushButton {"
692 "font-size:16px;"
693 "font-weight:bold;"
694 "border: 1px solid #000000;"
695 "border-radius: 9px;"
696 "color: %4;"
697 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
698 "}"
699 "QPushButton:hover {"
700 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
701 "border: 0px"
702 "}"
703 "QPushButton:flat {"
704 " border: none; /* no border for a flat push button */"
705 "}"
706 )
707 .arg( chgl->fopc.name())
708 .arg( chgl->fopc.darker(200).name())
709 .arg( chgl->fopc.lighter(100).name())
710 .arg("#ffffff"));
711
712 foMinusButton->setStyleSheet(QString(
713 "QPushButton {"
714 "font-size:16px;"
715 "font-weight:bold;"
716 "border: 1px solid #000000;"
717 "border-radius: 9px;"
718 "color: %4;"
719 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
720 "}"
721 "QPushButton:hover {"
722 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
723 "border: 0px"
724 "}"
725 "QPushButton:flat {"
726 " border: none; /* no border for a flat push button */"
727 "}"
728 )
729 .arg( chgl->fomc.name())
730 .arg( chgl->fomc.darker(200).name())
731 .arg( chgl->fomc.lighter(100).name())
732 .arg("#ffffff"));
733
734 diffMinusButton->setStyleSheet(QString(
735 "QPushButton {"
736 "font-size:16px;"
737 "font-weight:bold;"
738 "border: 1px solid #000000;"
739 "border-radius: 9px;"
740 "color: %4;"
741 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
742 "}"
743 "QPushButton:hover {"
744 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
745 "border: 0px"
746 "}"
747 "QPushButton:flat {"
748 " border: none; /* no border for a flat push button */"
749 "}"
750 )
751 .arg( chgl->dimc.name())
752 .arg( chgl->dimc.darker(200).name())
753 .arg( chgl->dimc.lighter(100).name())
754 .arg("#ffffff"));
755
756 diffPlusButton->setStyleSheet(QString(
757 "QPushButton {"
758 "font-size:16px;"
759 "font-weight:bold;"
760 "border: 1px solid #000000;"
761 "border-radius: 9px;"
762 "color: %4;"
763 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %3, stop: 1 %2);"
764 "}"
765 "QPushButton:hover {"
766 "background-color: qlineargradient(x1: 0, y1: 0, x2: 0, y2: 1, stop: 0 %1, stop: 0.5 %2, stop: 1 %1);"
767 "border: 0px"
768 "}"
769 "QPushButton:flat {"
770 " border: none; /* no border for a flat push button */"
771 "}"
772 )
773 .arg( chgl->dipc.name())
774 .arg( chgl->dipc.darker(200).name())
775 .arg( chgl->dipc.lighter(100).name())
776 .arg("#ffffff"));
777 mole->einstellung->beginGroup("Colors");
778 mole->einstellung->setValue("FourFoPlusColor", chgl->fopc);
779 mole->einstellung->setValue("FourFoMinusColor", chgl->fomc);
780 mole->einstellung->setValue("FourDiffPlusColor", chgl->dipc);
781 mole->einstellung->setValue("FourDiffMinusColor",chgl->dimc);
782 mole->einstellung->endGroup();
783 md->update();
784
785 }
786
loadFouAndPerform(const char filename[],bool neu)787 bool FourXle::loadFouAndPerform(const char filename[],bool neu){
788 /*! loads a fcf file prepares the reciprocal data for the fourier transpormation and performs it.
789 */
790 foti.start();
791
792 const int it[143]= {2,3,4,5,6,8,9,10,12,15,16,18,20,24,25,27,30,32,36,40,45,48,50,54,60,64,72,75,80,81,90,96,100,
793 108,120,125,128,135,144,150,160,162,180,192,200,216,225,240,243,250,256,270,288,300,320,324,360,375,384,400,405,
794 432,450,480,486,500,512,540,576,600,625,640,648,675,720,729,750,768,800,810,864,900,960,972,1000,1024,1080,1125,
795 1152,1200,1215,1250,1280,1296,1350,1440,1458,1500,1536,1600,1620,1728,1800,1875,1920,1944,2000,2025,2048,2160,
796 2187,2250,2304,2400,2430,2500,2560,2592,2700,2880,2916,3000,3072,3125,3200,3240,3375,3456,3600,3645,3750,3840,
797 3888,4000,4050,4096,4320,4374,4500,4608,4800,4860,5000};//!multiples of 2 3 and 5
798 mole->einstellung->beginGroup("Colors");
799 chgl->fopc=mole->einstellung->value("FourFoPlusColor", QColor("#0000FF")).value<QColor>();
800 chgl->fomc=mole->einstellung->value("FourFoMinusColor", QColor("#FFA600")).value<QColor>();
801 chgl->dipc=mole->einstellung->value("FourDiffPlusColor", QColor("#00B200")).value<QColor>();
802 chgl->dimc=mole->einstellung->value("FourDiffMinusColor",QColor("#CC0000")).value<QColor>();
803 mole->einstellung->endGroup();
804 //printf("%d\n",__LINE__);
805 killmaps();
806 //printf("%d {%s}\n",__LINE__,filename);
807 int ok;
808 if (!doMaps->isChecked()) return false;
809 if (strstr(filename,".fcf")==nullptr) return false;
810 FILE *f;
811 //printf("%d\n",__LINE__);
812 ok= readHeader(filename);
813 //printf("%d\n",__LINE__);
814 if (ok) {
815 switch (ok){
816 case 1: emit bigmessage("<font color=red>Map generation failed. SHELXL LIST code was not 6.</font>");break;
817 case 2: emit bigmessage(QString("<font color=red>Map generation failed. File %1 corrupted.</font>").arg(QString::fromLocal8Bit(filename)));break;
818 case 3: emit bigmessage(QString("<font color=red>Map generation failed. Cannot open file %1.</font>").arg(QString::fromLocal8Bit(filename)));break;
819 case 4: emit bigmessage(QString("<font color=red>Map generation failed. No reflection data in file %1.</font>").arg(QString::fromLocal8Bit(filename)));break;
820 }
821 return false;
822 }
823 //printf("%d\n",__LINE__);
824 f=fopen(filename,"rb");
825 if (f==nullptr)return false;
826 char line[122]="";
827 while (strstr(line,"_refln_phase_calc")==nullptr) {
828 if (fgets(line,120,f)){;}//just a stupid supression for annoying warnig
829 }
830 nr=0;
831 lr[nr].ih=0;
832 lr[nr].ik=0;
833 lr[nr].il=0;
834 lr[nr].fo=f000*f000;
835 lr[nr].so=0.5f;
836 lr[nr].fc=f000;
837 lr[nr].phi=0.0f;
838 nr=1;
839 //printf("f000 %g\n",f000);
840 emit bigmessage(QString::fromLocal8Bit(filename));
841
842 while (!feof(f)){
843 if (fgets(line,120,f)){;}//just a stupid supression for annoying warnig
844 int rdi=
845 sscanf(line,"%d %d %d %f %f %f %f",&lr[nr].ih,&lr[nr].ik, &lr[nr].il ,&lr[nr].fo, &lr[nr].so, &lr[nr].fc, &lr[nr].phi);
846 if (rdi==7) {
847 if ((abs(lr[nr].ih)<HKLMX)&&
848 (abs(lr[nr].ik)<HKLMX)&&
849 (abs(lr[nr].il)<HKLMX)&&
850 ((lr[nr].ih|lr[nr].ik|lr[nr].il)!=0))
851
852 //printf("%4d%4d%4d %g %g\n", lr[nr].ih, lr[nr].ik, lr[nr].il ,lr[nr].fo, lr[nr].so);
853 nr++;
854 }
855
856 // else printf("?? %d {%s}\n",rdi,line);
857
858 if (nr>=LM) {
859 emit bigmessage(QString("<font color=red>to many reflections in fcf file</font>"));
860 return false;
861 }
862 }
863 fclose(f);
864 if (nr<2) {
865 emit bigmessage(QString("<font color=red>Map generation failed. No reflection data in file %1.</font>").arg(QString::fromLocal8Bit(filename)));
866 return false;
867 }
868
869 for (int i=0;i<nr;i++){
870 double u=lr[i].ih,v=lr[i].ik,w=lr[i].il;
871 int mh=lr[i].ih,mk=lr[i].ik,ml=lr[i].il;
872 double p,q=lr[i].phi/180.0*M_PI;
873 lr[i].phi=fmod(4*M_PI+q,2*M_PI);
874 for (int k=0; k<ns; k++){
875 int nh,nk,nl;
876 double t=1.0;
877 nh=(int) (u*sy[0][k]+ v*sy[3][k] + w*sy[6][k]);
878 nk=(int) (u*sy[1][k]+ v*sy[4][k] + w*sy[7][k]);
879 nl=(int) (u*sy[2][k]+ v*sy[5][k] + w*sy[8][k]);
880 if((nl<0)||((nl==0)&&(nk<0))||((nl==0)&&(nk==0)&&(nh<0)))
881 {nh*=-1;nk*=-1;nl*=-1;t=-1.0;}
882 if ((nl<ml)||((nl==ml)&&(nk<mk))||((nl==ml)&&(nk==mk)&&(nh<=mh))) continue;
883 mh=nh;mk=nk;ml=nl;
884 p=u*sy[9][k]+v*sy[10][k]+w*sy[11][k];
885 lr[i].phi=fmod(4*M_PI+t*fmod(q-2*M_PI*p,2*M_PI)-0.01,2*M_PI)+0.01;
886
887 }
888 lr[i].ih=mh;
889 lr[i].ik=mk;
890 lr[i].il=ml;
891 }
892 sorthkl(nr,lr);
893 int n=-1;
894 {int i=0;
895 while(i<nr){
896 double t=0.;
897 double u=0.;
898 double v=0.;
899 double z=0.;
900 double p=0.;
901 int m;
902 int k=i;
903 while ((i<nr)&&(lr[i].ih==lr[k].ih)&&(lr[i].ik==lr[k].ik)&&(lr[i].il==lr[k].il)) {
904 t=t+1.;
905 u+=lr[i].fo;
906 v+=1./(lr[i].so*lr[i].so);
907 z+=lr[i].fc;
908 p=lr[i].phi;
909 i++;
910 }
911 m=n+1;
912 lr[m].fo=sqrt(fmax(0.,u/t));
913 lr[m].so=sqrt(lr[m].fo*lr[m].fo+sqrt(1./v))-lr[m].fo;
914 lr[m].fc=z/t;
915 lr[m].phi=p;
916 n=m;
917 lr[n].ih=lr[k].ih;
918 lr[n].ik=lr[k].ik;
919 lr[n].il=lr[k].il;
920 }
921 }
922 n++;
923 nr=n;
924 // printf("%4d%4d%4d %g %g %g %g\n",lr[0].ih,lr[0].ik,lr[0].il,lr[0].fo,lr[0].so,lr[0].fc,lr[0].phi);
925 {
926 float DX;
927 float DY;
928 float DZ;
929
930
931 {
932 int mh=0, mk=0, ml=0,j;
933 for (int n=0; n<nr; n++){
934 double u=lr[n].ih,v=lr[n].ik,w=lr[n].il;
935 int a,b,c;
936 for (int k=0; k<ns;k++){
937 a=abs((int)(u*sy[0][k]+v*sy[3][k]+w*sy[6][k]));
938 b=abs((int)(u*sy[1][k]+v*sy[4][k]+w*sy[7][k]));
939 c=abs((int)(u*sy[2][k]+v*sy[5][k]+w*sy[8][k]));
940 mh=(mh<a)?a:mh;
941 mk=(mk<b)?b:mk;
942 ml=(ml<c)?c:ml;
943 }
944 }
945 j=(int)(rr*mh+.5);
946 for (int i=0; it[i]< j; i++)n1=it[i+1];
947 j=(int)(rr*mk+.5);
948 for (int i=0; it[i]< j; i++)n2=it[i+1];
949 j=(int)(rr*ml+.5);
950 for (int i=0; (it[i]< j)||((nc)&&(it[i]%2)); i++) n3=it[i+1];
951 if (!voxelstr.isEmpty()) {
952 if (voxelstr.contains('x')) {
953 QStringList lll=voxelstr.split('x');
954 if (lll.size()==3){
955 n1=lll.at(0).toInt();
956 n2=lll.at(1).toInt();
957 n3=lll.at(2).toInt();
958 //printf("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n %d %d %d \n",n1,n2,n3);
959 }
960 }
961 }
962 n4=n2*n1;
963 n5=n3*n4;
964 datfo=(float*) malloc(sizeof(float)*n5);
965 if (datfo==nullptr){
966 fprintf(stderr,"Too less memory for F_observed!\n");
967 datfo=nullptr;
968 mapDefault();
969 info="";
970 emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
971 return false;
972 }
973 datfo_fc=(float*) malloc(sizeof(float)*n5);
974 if (datfo_fc==nullptr){
975 fprintf(stderr,"Too less memory for F_observed - F_calculated!\n");
976 free(datfo);
977 datfo=nullptr;
978 datfo_fc=nullptr;
979 mapDefault();
980 info="";
981 emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
982 return false;
983 }
984 DX=1.0f/n1;
985 DY=1.0f/n2;
986 DZ=1.0f/n3;
987 }
988 for (int typ=0; typ<2;typ++){
989 double miZ=99999.99,maZ=-99999.99;
990 // mInimum[typ]=99999.99; mAximum[typ]=-99999.99;
991 int nbytes,dims[3];
992 dims[0]=n3;
993 dims[1]=n2;
994 dims[2]=n1;
995 #ifdef FFTW3_H
996 B=(fftwf_complex*)fftwf_malloc(sizeof(fftwf_complex)*n5);
997 if (B==nullptr)return false;
998 for (int i=0; i<n5; i++){B[i][0]=0;B[i][1]=0;}
999 #else
1000 B=(kiss_fft_cpx*)KISS_FFT_MALLOC(nbytes = (sizeof(kiss_fft_cpx)*n5));
1001 if (B==nullptr){
1002 fprintf(stderr,"Too less memory for FFT map!\n");
1003 free(datfo);
1004 free(datfo_fc);
1005 datfo=nullptr;
1006 datfo_fc=nullptr;
1007 mapDefault();
1008 info="";
1009 emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
1010 return false;
1011 }
1012 for (int i=0; i<n5; i++){B[i].r=0;B[i].i=0;}
1013 #endif
1014 for (int i=0; i<nr;i++){
1015 float u,v,w;
1016 u=lr[i].ih;
1017 v=lr[i].ik;
1018 w=lr[i].il;
1019 float ss,s=0,t=0,q,p;
1020 for (int n=0; n<ns;n++){
1021 int j,k,l;
1022 j=(int) (u*sy[0][n]+ v*sy[3][n] + w*sy[6][n]);
1023 k=(int) (u*sy[1][n]+ v*sy[4][n] + w*sy[7][n]);
1024 l=(int) (u*sy[2][n]+ v*sy[5][n] + w*sy[8][n]);
1025 if((abs(j-lr[i].ih)+abs(k-lr[i].ik)+abs(l-lr[i].il))==0)s+=1.0;
1026 if(abs(j+lr[i].ih)+abs(k+lr[i].ik)+abs(l+lr[i].il)==0)t+=1.0;
1027 }
1028 if (i==0) {//printf("v%f s%f t%f\n",C[14],s,t);
1029 s=1;t=0;//f000
1030 }
1031 if(typ==0) ss=(lr[i].fo-lr[i].fc)/(C[14]*(s+t));
1032 else if(typ==2) //here this happens never
1033 ss=(lr[i].fc)/(C[14]*(s+t));
1034 else ss=(lr[i].fo)/(C[14]*(s+t));
1035 if(lr[i].fc>1.E-6) ss=ss/(1.+rw*pow(lr[i].so/lr[i].fc,4));
1036 for (int n=0; n<ns;n++){
1037 int j,k,l,m;
1038 j=(int) (u*sy[0][n]+ v*sy[3][n] + w*sy[6][n]);
1039 k=(int) (u*sy[1][n]+ v*sy[4][n] + w*sy[7][n]);
1040 l=(int) (u*sy[2][n]+ v*sy[5][n] + w*sy[8][n]);
1041 // q=(-2*M_PI*(u*sy[9][n]+v*sy[10][n]+w*sy[11][n]))-M_PI*(j*DX+k*DY+l*DZ);
1042 q=(lr[i].phi-2*M_PI*(u*sy[9][n]+v*sy[10][n]+w*sy[11][n]))-M_PI*(j*DX+k*DY+l*DZ);
1043 j=(999*n1+j)%n1;
1044 k=(999*n2+k)%n2;
1045 l=(999*n3+l)%n3;
1046 m=j+n1*(k+n2*l);
1047 p=ss*cosf(q);
1048 q=ss*sinf(q);
1049 #ifdef FFTW3_H
1050 B[m][0]=p;
1051 B[m][1]=q;
1052 #else
1053 B[m].r=p;
1054 B[m].i=q;
1055 #endif
1056 j*=-1;
1057 if(j<0)j=n1+j;
1058 k*=-1;
1059 if(k<0)k=n2+k;
1060 l*=-1;
1061 if(l<0)l=n3+l;
1062 m=j+n1*(k+n2*l);
1063 #ifdef FFTW3_H
1064 B[m][0]=p;
1065 B[m][1]=-q;
1066 #else
1067 B[m].r=p;
1068 B[m].i=-q;
1069 #endif
1070 }
1071 }
1072 #ifdef FFTW3_H
1073 fwd_plan = fftwf_plan_dft_3d(n3,n2,n1,B,B,FFTW_FORWARD,FFTW_ESTIMATE);
1074 fftwf_execute(fwd_plan);
1075 fftwf_destroy_plan(fwd_plan);
1076 #else
1077 fwd_plan = kiss_fftnd_alloc(dims,3,0,0,0);
1078 kiss_fftnd( fwd_plan,B,B);
1079 free(fwd_plan);
1080 #endif
1081 float t=0;
1082 double DM=0., DS=0., DD ;
1083 for (int i=0; i<n5;i++){
1084 #ifdef FFTW3_H
1085 DD=B[i][0];
1086 #else
1087 DD=B[i].r;
1088 #endif
1089 // maxi[typ]=fmax(maxi[typ],DD); mini[typ]=fmin(mini[typ],DD);
1090 miZ=fmin(miZ,DD);
1091 maZ=fmax(maZ,DD);
1092 DM+=DD;
1093 DS+=DD*DD;
1094 #ifdef FFTW3_H
1095 if (typ==1) datfo[i]=B[i][0];
1096 else if (typ==0) datfo_fc[i]=B[i][0];
1097 #else
1098 if (typ==1) datfo[i]=B[i].r;
1099 else if (typ==0) datfo_fc[i]=B[i].r;
1100
1101 #endif
1102 }
1103 sigma[typ]=t=sqrt((DS/n5)-((DM/n5)*(DM/n5)));
1104 free(B);
1105 }//1
1106 }//2
1107 urs=V3(0,0,0);int gt=0;
1108 for (int i=0; i<mole->showatoms.size();i++)
1109 if (mole->showatoms.at(i).an>-1) {urs+=mole->showatoms.at(i).frac;gt++;}
1110 urs*=1.0/gt;
1111 urs=V3(1,1,1)-1.0*urs;
1112 mole->frac2kart(urs,urs);
1113
1114 nodex= (FNode*)malloc(sizeof(FNode)*n5);
1115 if (nodex==nullptr){
1116 fprintf(stderr,"Too less memory for node X!\n");
1117 free(datfo);
1118 free(datfo_fc);
1119 datfo=nullptr;
1120 datfo_fc=nullptr;
1121 mapDefault();
1122 info="";
1123 emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
1124 return false;
1125 }
1126 nodey= (FNode*)malloc(sizeof(FNode)*n5);
1127 if (nodey==nullptr){
1128 fprintf(stderr,"Too less memory for node Y!\n");
1129 free(nodex);
1130 free(datfo);
1131 free(datfo_fc);
1132 nodex = nullptr;
1133 datfo = nullptr;
1134 datfo_fc = nullptr;
1135 mapDefault();
1136 info="";
1137 emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
1138 return false;
1139 }
1140 nodez= (FNode*)malloc(sizeof(FNode)*n5);
1141 if (nodez==nullptr){
1142 fprintf(stderr,"Too less memory for node Z!\n");
1143 free(nodex);
1144 free(nodey);
1145 free(datfo);
1146 free(datfo_fc);
1147 nodex = nullptr;
1148 nodey = nullptr;
1149 datfo = nullptr;
1150 datfo_fc = nullptr;
1151 mapDefault();
1152 info="";
1153 emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
1154 return false;
1155 }
1156 noGoMap= (int*) malloc(sizeof(int)*n5*27);
1157 if (noGoMap==nullptr){
1158 fprintf(stderr,"Too less memory for noGoMap!\n");
1159
1160 free(nodex);
1161 free(nodey);
1162 free(nodez);
1163 free(datfo);
1164 free(datfo_fc);
1165 nodex= nullptr;
1166 nodey= nullptr;
1167 nodez= nullptr;
1168 datfo= nullptr;
1169 datfo_fc= nullptr;
1170 mapDefault();
1171 info="";
1172 emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
1173 return false;
1174 }
1175 for (int no=0;no<(n5*27);no++) noGoMap[no]=0;
1176 for (int o=0; o<n5;o++){
1177 nodex[o].flag=0;
1178 nodey[o].flag=0;
1179 nodez[o].flag=0;
1180 }
1181 dx=V3(1.0/(n1),0,0);
1182 dy=V3(0,1.0/(n2),0);
1183 dz=V3(0,0,1.0/(n3));
1184 mole->frac2kart(dx,dx);
1185 mole->frac2kart(dy,dy);
1186 mole->frac2kart(dz,dz);
1187 delDA[ 0]= -n1*dx -n2*dy -n3*dz;//nx ny,nz??
1188 delDA[ 1]= -n2*dy -n3*dz;
1189 delDA[ 2]= n1*dx -n2*dy -n3*dz;
1190 delDA[ 3]= -n1*dx -n3*dz;
1191 delDA[ 4]= -n3*dz;
1192 delDA[ 5]= n1*dx -n3*dz;
1193 delDA[ 6]= -n1*dx +n2*dy -n3*dz;
1194 delDA[ 7]= n2*dy -n3*dz;
1195 delDA[ 8]= n1*dx +n2*dy -n3*dz;
1196 delDA[ 9]= -n1*dx -n2*dy ;
1197 delDA[10]= -n2*dy ;
1198 delDA[11]= n1*dx -n2*dy ;
1199 delDA[12]= -n1*dx ;
1200 delDA[13]= V3(0,0,0) ;
1201 delDA[14]= n1*dx ;
1202 delDA[15]= -n1*dx +n2*dy ;
1203 delDA[16]= n2*dy ;
1204 delDA[17]= n1*dx +n2*dy ;
1205 delDA[18]= -n1*dx -n2*dy +n3*dz;
1206 delDA[19]= -n2*dy +n3*dz;
1207 delDA[20]= n1*dx -n2*dy +n3*dz;
1208 delDA[21]= -n1*dx +n3*dz;
1209 delDA[22]= +n3*dz;
1210 delDA[23]= n1*dx +n3*dz;
1211 delDA[24]= -n1*dx +n2*dy +n3*dz;
1212 delDA[25]= n2*dy +n3*dz;
1213 delDA[26]= n1*dx +n2*dy +n3*dz;
1214 emit bigmessage(QString("Uniq Reflections: %1 <br>Fourier grid dimensions: %2 %3 %4 = %5 grid points.<br>Time for loading and fourier transformation: <b>%6 s</b>.<br>").arg(nr).arg(n1).arg(n2).arg(n3).arg(n5).arg(foti.elapsed()/1000.0,1,'f',1));
1215 if (neu) gen_surface(neu);
1216 return true;
1217 }
1218
loadMap(int NX,int NY,int NZ,float * dat)1219 bool FourXle::loadMap(int NX, int NY, int NZ, float *dat){
1220 bool wbi=isBelo;
1221 //printf("%d\n",__LINE__);
1222 killmaps();
1223 //printf("%d\n",__LINE__);
1224 isBelo=wbi;
1225
1226 n1=NX;
1227 n2=NY;
1228 n3=NZ;
1229 n4=n2*n1;
1230 n5=n3*n4;
1231 datfo=(float*) malloc(sizeof(float)*n5);
1232 if (datfo==nullptr){
1233 fprintf(stderr,"Too less memory for F_observed - F_calculated!\n");
1234 free(datfo_fc);
1235 mapDefault();
1236 info="";
1237 emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
1238 return false;
1239 }
1240 for (int i=0; i<n5; i++)datfo[i]=0.0f;
1241 datfo_fc=dat;
1242
1243 sigma[0] = sigmaFoFc();
1244 sigma[1] = sigmaFo();
1245 iso[1]= sigma[0] * 2.5f;
1246 iso[0]= sigma[1] * 1.2f;
1247 fomaps->setValue(iso[0]);
1248 difmaps->setValue(iso[1]);
1249
1250 urs=V3(0,0,0);int gt=0;
1251 for (int i=0; i<mole->showatoms.size();i++)
1252 if (mole->showatoms.at(i).an>-1) {urs+=mole->showatoms.at(i).frac;gt++;}
1253 urs*=1.0/gt;
1254 urs=V3(1,1,1)-1.0*urs;
1255 mole->frac2kart(urs,urs);
1256 nodex= (FNode*)malloc(sizeof(FNode)*n5);
1257 if (nodex==nullptr){
1258 fprintf(stderr,"Too less memory for node X!\n");
1259 free(datfo);
1260 free(datfo_fc);
1261 mapDefault();
1262 info="";
1263 emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
1264 return false;
1265 }
1266 nodey= (FNode*)malloc(sizeof(FNode)*n5);
1267 if (nodey==nullptr){
1268 fprintf(stderr,"Too less memory for node Y!\n");
1269 free(datfo);
1270 free(datfo_fc);
1271 free(nodex);
1272 mapDefault();
1273 info="";
1274 emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
1275 return false;
1276 }
1277 nodez= (FNode*)malloc(sizeof(FNode)*n5);
1278 if (nodez==nullptr){
1279 fprintf(stderr,"Too less memory for node Z!\n");
1280 free(datfo);
1281 free(datfo_fc);
1282 free(nodex);
1283 free(nodey);
1284 mapDefault();
1285 info="";
1286 emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
1287 return false;
1288 }
1289 noGoMap= (int*) malloc(sizeof(int)*n5*27);
1290 if (noGoMap==nullptr){
1291 fprintf(stderr,"Too less memory for noGoMap!\n");
1292 free(datfo);
1293 free(datfo_fc);
1294 free(nodex);
1295 free(nodey);
1296 free(nodez);
1297 mapDefault();
1298 info="";
1299 emit bigmessage(QString("<font color=red>Could not allocate sufficient memory for density maps. Defaults have been applied.</font> (Try to tick and un tick fcf6 by hkl)"));
1300 return false;
1301 }
1302 for (int no=0;no<(n5*27);no++) noGoMap[no]=0;
1303 for (int o=0; o<n5;o++){
1304 nodex[o].flag=0;
1305 nodey[o].flag=0;
1306 nodez[o].flag=0;
1307 }
1308 dx = V3(1.0 / (n1), 0, 0);
1309 dy = V3(0, 1.0 / (n2), 0);
1310 dz = V3(0, 0, 1.0 / (n3));
1311 mole->frac2kart(dx,dx);
1312 mole->frac2kart(dy,dy);
1313 mole->frac2kart(dz,dz);
1314 delDA[ 0]= -n1*dx -n2*dy -n3*dz;//nx ny,nz??
1315 delDA[ 1]= -n2*dy -n3*dz;
1316 delDA[ 2]= n1*dx -n2*dy -n3*dz;
1317 delDA[ 3]= -n1*dx -n3*dz;
1318 delDA[ 4]= -n3*dz;
1319 delDA[ 5]= n1*dx -n3*dz;
1320 delDA[ 6]= -n1*dx +n2*dy -n3*dz;
1321 delDA[ 7]= n2*dy -n3*dz;
1322 delDA[ 8]= n1*dx +n2*dy -n3*dz;
1323 delDA[ 9]= -n1*dx -n2*dy ;
1324 delDA[10]= -n2*dy ;
1325 delDA[11]= n1*dx -n2*dy ;
1326 delDA[12]= -n1*dx ;
1327 delDA[13]= V3(0,0,0) ;
1328 delDA[14]= n1*dx ;
1329 delDA[15]= -n1*dx +n2*dy ;
1330 delDA[16]= n2*dy ;
1331 delDA[17]= n1*dx +n2*dy ;
1332 delDA[18]= -n1*dx -n2*dy +n3*dz;
1333 delDA[19]= -n2*dy +n3*dz;
1334 delDA[20]= n1*dx -n2*dy +n3*dz;
1335 delDA[21]= -n1*dx +n3*dz;
1336 delDA[22]= +n3*dz;
1337 delDA[23]= n1*dx +n3*dz;
1338 delDA[24]= -n1*dx +n2*dy +n3*dz;
1339 delDA[25]= n2*dy +n3*dz;
1340 delDA[26]= n1*dx +n2*dy +n3*dz;
1341 rr = mapprec->value();
1342 rw = weak->value();
1343 //gen_surface_silent();
1344 return true;
1345
1346 }
1347
1348
trimm(char s[])1349 void FourXle::trimm(char s[]){
1350 /*! a trimm function for c-strings.
1351 */
1352 char sc[409];
1353 int j=0;
1354 int len=strlen(s);
1355 strncpy(sc,s,400);
1356 for (int i=0; i<len; i++) if ((sc[i]!='\'')&&(!isspace(sc[i]))) s[j++]=toupper(sc[i]);
1357 s[j]='\0';
1358 }
1359
deletes(char * s,int count)1360 void FourXle::deletes(char *s, int count){
1361 /*! deletes count characters at the begining of s.
1362 */
1363 if ((s==nullptr)||(count <1)||((size_t)count>strlen(s))) return;
1364 for (int i=0; i<count; i++) s[i]=' ';
1365 trimm(s);
1366 }
1367
exportMaps(int na,const char filename[],const char atomlist[])1368 void FourXle::exportMaps(int na, const char filename[], const char atomlist[]){
1369 if ((datfo==nullptr)||(datfo_fc==nullptr)) return;
1370 FILE *fo,*fof1;//,*f1f2;
1371 char foname[4096];
1372 char fof1name[4096];
1373 // char f1f2name[4096];
1374 int len=strlen(filename);
1375 float factor=0.1481847095290449f;//a0**3
1376 double a0=0.52917720859; //bohr=0.5291775108
1377 int i=0;
1378 //
1379 //FO MAP
1380 //
1381 V3 dx1=V3(1.0/(n1),0,0);
1382 V3 dy1=V3(0,1.0/(n2),0);
1383 V3 dz1=V3(0,0,1.0/(n3));
1384 mole->frac2kart(dx1,dx1);
1385 mole->frac2kart(dy1,dy1);
1386 mole->frac2kart(dz1,dz1);
1387
1388 strncpy(foname,filename,len-4);
1389 foname[len-4]='\0';
1390 strcat(foname,"_fo_densitymap.cube");
1391 fo=fopen(foname,"w");
1392 if (fo==nullptr) return ;
1393 double fomax=-9e37,fomin=9e37,em=0,ep=0,et=0;
1394 for (int i=0; i<n5; i++ ){
1395 if (datfo[i]>0.0) ep+=datfo[i];
1396 else em+=datfo[i];
1397 et+=datfo[i];
1398 fomax=fmax(datfo[i],fomax);
1399 fomin=fmin(datfo[i],fomin);
1400 }
1401 //printf("Fo %g %g %g %g %g\n",fomin,fomax,ep/n5*C[14],em/n5*C[14],et/n5*C[14] );
1402 fprintf(fo,"F observed map written by ShelXle\nDensity obtained from Fo in fcf with phases of model transformed using fft\n");
1403 fprintf(fo,"%5d%12.6f%12.6f%12.6f\n",na,
1404 (dx1.x+dy1.x+dz1.x)*0.5/a0,
1405 (dx1.y+dy1.y+dz1.y)*0.5/a0,
1406 (dx1.z+dy1.z+dz1.z)*0.5/a0
1407 //0.0,0.0,0.0
1408 );
1409 fprintf(fo,"%5d%12.6f%12.6f%12.6f\n",n1,dx1.x/a0,dx1.y/a0,dx1.z/a0);
1410 fprintf(fo,"%5d%12.6f%12.6f%12.6f\n",n2,dy1.x/a0,dy1.y/a0,dy1.z/a0);
1411 fprintf(fo,"%5d%12.6f%12.6f%12.6f" ,n3,dz1.x/a0,dz1.y/a0,dz1.z/a0);//no newline here because it is in atomlist
1412 fprintf(fo,"%s",atomlist);
1413 for (int xi=0;xi<n1;xi++)
1414 for (int yi=0;yi<n2;yi++)
1415 for (int zi=0;zi<n3;zi++)
1416 {
1417 //fprintf(fo,"%s%13.5E",(((i%6)==0)?"\n":""),datfo[dex(xi,yi,zi)]*factor);
1418 fprintf(fo,"%s",QString("%1%2").arg(((i%6)==0)?"\n":"").arg(datfo[dex(xi,yi,zi)]*factor,13,'E',5).toStdString().c_str());
1419 i++;
1420 }
1421 fclose(fo);
1422 //
1423 //FO-F1 MAP
1424 //
1425
1426 strncpy(fof1name,filename,len-4);
1427 fof1name[len-4]='\0';
1428 strcat(fof1name,"_fo-fc_densitymap.cube");
1429 fof1=fopen(fof1name,"w");
1430 if (fof1==nullptr) return ;
1431 fomax=-9e37;
1432 fomin= 9e37;
1433 em=0;
1434 ep=0;
1435 et=0;
1436 for (int i=0; i<n5; i++ ){
1437 if (datfo_fc[i]>0.0) ep+=datfo_fc[i];
1438 else em+=datfo_fc[i];
1439 et+=datfo_fc[i];
1440 fomax=fmax(datfo_fc[i],fomax);
1441 fomin=fmin(datfo_fc[i],fomin);
1442 }
1443 printf("Fo-Fc %g %g %g %g %g\n",fomin,fomax,ep/n5*C[14],em/n5*C[14],et/n5*C[14] );
1444 fprintf(fof1,"F obs-fc map written by ShelXle\nDifferene density Fo-fc in fcf file transformed using fft\n");
1445 fprintf(fof1,"%5d%12.6f%12.6f%12.6f\n",na,
1446 (dx1.x+dy1.x+dz1.x)*0.5/a0,
1447 (dx1.y+dy1.y+dz1.y)*0.5/a0,
1448 (dx1.z+dy1.z+dz1.z)*0.5/a0
1449 //0.0,0.0,(dx.z+dy.z+dz.z)/a0
1450 );
1451 fprintf(fof1,"%5d%12.6f%12.6f%12.6f\n",n1,dx1.x/a0,dx1.y/a0,dx1.z/a0);
1452 fprintf(fof1,"%5d%12.6f%12.6f%12.6f\n",n2,dy1.x/a0,dy1.y/a0,dy1.z/a0);
1453 fprintf(fof1,"%5d%12.6f%12.6f%12.6f" ,n3,dz1.x/a0,dz1.y/a0,dz1.z/a0);//no newline here because it is in atomlist
1454 fprintf(fof1,"%s",atomlist);
1455 for (int xi=0;xi<n1;xi++)
1456 for (int yi=0;yi<n2;yi++)
1457 for (int zi=0;zi<n3;zi++)
1458 {
1459 //fprintf(fof1,"%s%13.5E",(((i%6)==0)?"\n":""),datfo_fc[dex(xi,yi,zi)]*factor);
1460 fprintf(fof1,"%s",QString("%1%2").arg(((i%6)==0)?"\n":"").arg(datfo_fc[dex(xi,yi,zi)]*factor,13,'E',5).toStdString().c_str());
1461 i++;
1462 }
1463 fclose(fof1);
1464 }
1465 // */
1466 //}
1467
readHeader(const char * filename)1468 int FourXle::readHeader(const char *filename){
1469 /*! reads the header of an fcf file
1470 */
1471 FILE *f=nullptr;
1472 char line[122],*dum=nullptr;
1473 //size_t zlen=120;
1474 int ok=0;
1475 int i;double T,V;
1476 f=fopen(filename,"r");
1477 //printf("%d\n",__LINE__);
1478 if (f==nullptr) return 3;
1479 //printf("%d\n",__LINE__);
1480 ns=0;
1481 sy[0][ns]=1.0;
1482 sy[1][ns]=0.0;
1483 sy[2][ns]=0.0;
1484
1485 sy[3][ns]=0.0;
1486 sy[4][ns]=1.0;
1487 sy[5][ns]=0.0;
1488
1489 sy[6][ns]=0.0;
1490 sy[7][ns]=0.0;
1491 sy[8][ns]=1.0;
1492
1493 sy[9][ns]=0.0;
1494 sy[10][ns]=0.0;
1495 sy[11][ns]=0.0;
1496 ns=1;
1497 int listcode=0;
1498 do{
1499 dum=fgets(line,120,f);
1500 //printf("%d\n",__LINE__);
1501 if (dum==nullptr){fclose(f);return 2;};
1502 if (feof(f)){fclose(f);return 2;};
1503 while (dum[0]==' ') dum++;
1504 if (!strncmp(dum,"_shelx_title",12)) {
1505 sscanf(line,"_shelx_title %[^\r\n]",titl);
1506 trimm(titl);
1507 }
1508 if (!strncmp(dum,"_exptl_crystal_F_000",20)){
1509 sscanf(line,"_exptl_crystal_F_000 %f",&f000);
1510 if (overridef000!=-1.0) f000=overridef000;
1511 //printf("F000 = %g %g %d %s\n",f000,overridef000,(overridef000!=-1.0f),line);
1512 }
1513 if (!strncmp(dum,"_shelx_refln_list_code",22)) {
1514 sscanf(line,"_shelx_refln_list_code %d",&listcode);
1515 //qDebug()<<listcode;
1516 if (listcode!=6) {fclose(f);return 1;}
1517 }
1518 if (!strncmp(dum,"_cell_length_a",14)) {
1519 sscanf(line,"_cell_length_a %lf",&C[0]);
1520 }
1521 if (!strncmp(dum,"_cell_length_b",14)) {
1522 sscanf(line,"_cell_length_b %lf",&C[1]);
1523 }
1524 if (!strncmp(dum,"_cell_length_c",14)) {
1525 sscanf(line,"_cell_length_c %lf",&C[2]);
1526 }
1527 if (!strncmp(dum,"_cell_angle_alpha",17)) {
1528 sscanf(line,"_cell_angle_alpha %lf",&C[3]);
1529 }
1530 if (!strncmp(dum,"_cell_angle_beta",16)) {
1531 sscanf(line,"_cell_angle_beta %lf",&C[4]);
1532 }
1533 if (!strncmp(dum,"_cell_angle_gamma",17)) {
1534 sscanf(line,"_cell_angle_gamma %lf",&C[5]);
1535 for (i=0;i<3;i++){
1536 if (C[i]<0.1) return 2;
1537 T=.0174533*C[i+3];
1538 if (T<0.001) return 2;
1539 D[i]=sin(T);
1540 D[i+3]=cos(T);
1541 C[i+6]=(D[i]/(C[i]*C[i]));
1542 }
1543 V=1.-D[3]*D[3]-D[4]*D[4]-D[5]*D[5]+2.*D[3]*D[4]*D[5];
1544 C[6]/=V;
1545 C[7]/=V;
1546 C[8]/=V;
1547 C[9]= 2.*sqrt(C[7]*C[8])*(D[4]*D[5]-D[3])/(D[2]*D[2]);
1548 C[10]=2.*sqrt(C[6]*C[8])*(D[3]*D[5]-D[4])/(D[0]*D[2]);
1549 C[11]=2.*sqrt(C[6]*C[7])*(D[3]*D[4]-D[5])/(D[0]*D[1]);
1550 C[14]=C[1]*C[2]*C[0]*sqrt(V);
1551 D[6]=C[1]*C[2]*D[0]/C[14];
1552 D[7]=C[0]*C[2]*D[1]/C[14];
1553 D[8]=C[0]*C[1]*D[2]/C[14];
1554 }
1555 if ((!strncmp(dum,"_symmetry_equiv_pos_as_xyz",26))||(!strncmp(dum,"_space_group_symop_operation_xyz",32))){
1556 // char s1[50],s2[50],s3[50];
1557 // char *kill=nullptr,*nom=nullptr,*div=nullptr ;
1558 dum=fgets(line,120,f);
1559 trimm(line);
1560 while (strchr(line,'Y')) {
1561 QString sc=QString(line).toUpper().remove("SYMM").trimmed();
1562 sc=sc.remove("'");
1563 sc=sc.remove(" ");
1564 QStringList axe=sc.split(",");
1565 QStringList bruch;
1566 if (axe.size()!=3) return false;
1567 double _sx[3],_sy[3],_sz[3],t[3];
1568 for (int i=0; i<3; i++){
1569 _sx[i]=0;_sy[i]=0;_sz[i]=0;t[i]=0;
1570 if (axe.at(i).contains("-X")) {_sx[i]=-1.0;axe[i].remove("-X");}
1571 else if (axe.at(i).contains("X")) {_sx[i]=1.0;axe[i].remove("X");}
1572 if (axe.at(i).contains("-Y")) {_sy[i]=-1.0;axe[i].remove("-Y");}
1573 else if (axe.at(i).contains("Y")) {_sy[i]=1.0;axe[i].remove("Y");}
1574 if (axe.at(i).contains("-Z")) {_sz[i]=-1.0;axe[i].remove("-Z");}
1575 else if (axe.at(i).contains("Z")) {_sz[i]=1.0;axe[i].remove("Z");}
1576 if (axe.at(i).endsWith("+")) axe[i].remove("+");
1577 if (axe.at(i).contains("/")) {
1578 bruch=axe.at(i).split("/");
1579 if (bruch.size()==2) t[i]=bruch.at(0).toDouble() / bruch.at(1).toDouble();
1580 }
1581 else if (!axe.at(i).isEmpty()) t[i]=axe.at(i).toDouble();
1582 }
1583 sy[0][ns]=_sx[0];
1584 sy[1][ns]=_sy[0];
1585 sy[2][ns]=_sz[0];
1586
1587 sy[3][ns]=_sx[1];
1588 sy[4][ns]=_sy[1];
1589 sy[5][ns]=_sz[1];
1590
1591 sy[6][ns]=_sx[2];
1592 sy[7][ns]=_sy[2];
1593 sy[8][ns]=_sz[2];
1594
1595
1596 sy[9][ns]=t[0];
1597 sy[10][ns]=t[1];
1598 sy[11][ns]=t[2];
1599
1600
1601 strcpy(line,"");
1602 dum=fgets(line,120,f);
1603 trimm(line);
1604 ns++;
1605
1606 }
1607 }
1608 if (!strncmp(dum,"_refln_phase_calc",17)) ok=1;
1609 }while((!ok)&&(!feof(f)));
1610
1611 if (listcode!=6) return 1;
1612 for (int i=0; i<ns;i++){
1613 for (int n=i+1; n<ns;n++){
1614 int u=0,v=0;
1615 for (int j=0; j<9; j++){
1616 u+=abs(sy[j][n]-sy[j][i]);
1617 v+=abs(sy[j][n]+sy[j][i]);
1618 }
1619 if (fmin(u,v)>0.01) continue;
1620 for (int j=0; j<12; j++){
1621 sy[j][n]=sy[j][ns-1];
1622 }
1623 ns--;
1624 }
1625 }
1626 fclose(f);
1627
1628 return 0;
1629 }
1630
sorthkl(int nr,Rec r[])1631 void FourXle::sorthkl(int nr, Rec r[]){
1632 /*! sorts the reflection list
1633 */
1634 Rec *hilf= (Rec*) malloc(sizeof(Rec)*nr);
1635 if (hilf==nullptr)return ;
1636 int i,j,k,nj,ni,spalte;int index[4096];
1637 for (spalte=0; spalte<3; spalte++){
1638 j=-999999;
1639 k=999999;
1640 switch (spalte) {
1641 case 0: for (i=0; i<nr; i++){ j=(j<r[i].ih)?r[i].ih:j; k=(k>r[i].ih)?r[i].ih:k; } break;
1642 case 1: for (i=0; i<nr; i++){ j=(j<r[i].ik)?r[i].ik:j; k=(k>r[i].ik)?r[i].ik:k; } break;
1643 case 2: for (i=0; i<nr; i++){ j=(j<r[i].il)?r[i].il:j; k=(k>r[i].il)?r[i].il:k; } break;
1644 }
1645 nj=-k;
1646 ni=(nj+j+1);
1647 for (i=0; i<=ni; i++) index[i]=0;
1648 for (i=0; i<nr; i++){
1649 switch (spalte){
1650 case 0: j=r[i].ih+nj; break;
1651 case 1: j=r[i].ik+nj; break;
1652 case 2: j=r[i].il+nj; break;
1653 }
1654 index[j]++;/*brauch ich das? -->JA!*/
1655 hilf[i].ih=r[i].ih;
1656 hilf[i].ik=r[i].ik;
1657 hilf[i].il=r[i].il;
1658 hilf[i].fo=r[i].fo;
1659 hilf[i].so=r[i].so;
1660 hilf[i].fc=r[i].fc;
1661 hilf[i].phi=r[i].phi;
1662 }/*/4*/
1663 j=0;
1664 for (i=0; i<ni; i++){
1665 k=j;
1666 j+=index[i];
1667 index[i]=k;
1668 }/*/5*/
1669 for (i=0; i<nr;i++){
1670 switch (spalte) {
1671 case 0: j=hilf[i].ih +nj;break;
1672 case 1: j=hilf[i].ik +nj;break;
1673 case 2: j=hilf[i].il +nj;break;
1674 }
1675 index[j]++;
1676 j=index[j]-1;
1677 r[j].ih=hilf[i].ih;
1678 r[j].ik=hilf[i].ik;
1679 r[j].il=hilf[i].il;
1680 r[j].fo=hilf[i].fo;
1681 r[j].so=hilf[i].so;
1682 r[j].fc=hilf[i].fc;
1683 r[j].phi=hilf[i].phi;
1684 }/*/6*/
1685 }/*/spalten*/
1686 free(hilf);
1687 }
1688
bewegt(V3 nm)1689 void FourXle::bewegt(V3 nm){/*!moves the rotation center to nm*/
1690 // printf("moved\n");
1691 V3 v , alturs=urs;
1692 mole->kart2frac(nm,v);
1693 urs=V3(1,1,1)-1.0*v;
1694 mole->frac2kart(urs,urs);
1695 if ((chgl->objCnt==acnt)&&(maptrunc==2)) return;
1696 if (urs==alturs) return;
1697 gen_surface(false);
1698 chgl->pause=false;
1699 chgl->updateGL();
1700 }
1701
inimap()1702 void FourXle::inimap(){/*! reinitializes the display lists for screenshots*/
1703 gen_surface(false);
1704
1705 }
1706
gen_surface(bool neu,int imin,int imax)1707 void FourXle::gen_surface(bool neu,int imin,int imax){
1708 GenMapWorker *w = new GenMapWorker(chgl, mole, this, neu, imin, imax);
1709 chgl->habzutun=true;
1710 //connect(this,SIGNAL(haltMW()),w,SLOT(terminate()));
1711 connect(w, SIGNAL(finished()),this, SLOT(reportInfo()));
1712 connect(w, SIGNAL(finished()), chgl, SLOT(updateGL()));
1713 connect(w, SIGNAL(finished()), w, SLOT(deleteLater()));
1714 connect(w, SIGNAL(finished()),chgl, SLOT(fertig()));
1715 w->start();
1716 }
1717 //GenMapWorker
1718
1719
1720
gen_surface(bool neu,int imin,int imax)1721 void GenMapWorker::gen_surface(bool neu,int imin,int imax){
1722 if (fxle->noGoMap==nullptr) return;
1723 /*! creates iso surfaces for fo-fc- fo-fc+ fo+ and fo- maps if neu then the iso values are calculated fro the sigma value of the map.
1724 *
1725 */
1726 if (!chgl->neutrons) imax=qMin(3,imax);
1727 V3 v=V3(0,0,0);
1728 mole->kart2frac(chgl->altemitte,v);
1729 fxle->urs=V3(1,1,1)-1.0*v;
1730 mole->frac2kart(fxle->urs,fxle->urs);
1731 //chgl->pause=true;
1732 chgl->fVertexes.resize(28);
1733 chgl->fNormals.resize(28);
1734 fxle->suti.start();
1735 fxle->disconnect(chgl,SIGNAL(diffscroll(int ,int )),0,0);
1736 fxle->disconnect(chgl,SIGNAL(neuemitte(V3)),0,0);
1737 fxle->disconnect(chgl,SIGNAL(inimibas()),0,0);
1738 // if ((neu)||(mole->showatoms.size()!=oldatomsize)||(chgl->unhide->isVisible())){
1739 for (int no=0;no<(fxle->n5*27);no++) fxle->noGoMap[no]=0;
1740 mole->kart2frac(fxle->dx,fxle->dxc);
1741 mole->kart2frac(fxle->dy,fxle->dyc);
1742 mole->kart2frac(fxle->dz,fxle->dzc);
1743 int incx = (int) (fxle->cutrange->value()/fxle->dx.x);
1744 int incy = (int) (fxle->cutrange->value()/fxle->dy.y);
1745 int incz = (int) (fxle->cutrange->value()/fxle->dz.z);
1746 int incmin=qMin(incx,qMin(incy,incz));
1747 incmin*=incmin;
1748 for (int g=0; g<mole->showatoms.size();g++){
1749 if (mole->showatoms.at(g).hidden) continue;
1750 mole->kart2frac(mole->showatoms.at(g).pos,fxle->oc);
1751 int ax=(int)((fxle->oc.x)/fxle->dxc.x-0.499), bx=(int)((fxle->oc.y)/fxle->dyc.y-0.499), cx=(int)((fxle->oc.z)/fxle->dzc.z-0.499);
1752 for (int aa=ax-incx; aa<ax+incx; aa++){
1753 for (int bb=bx-incy; bb<bx+incy; bb++){
1754 for (int cc=cx-incz; cc<cx+incz; cc++){
1755 if (incmin<((aa-ax)*(aa-ax)+(bb-bx)*(bb-bx)+(cc-cx)*(cc-cx))) continue;
1756 int dEx=fxle->dex3(aa,bb,cc);
1757 if ((dEx>0)&&(dEx<(27*fxle->n5-1)))
1758 fxle->noGoMap[dEx]=1;
1759 }
1760 }
1761 }
1762 }
1763 // oldatomsize=mole->showatoms.size();
1764 // }
1765 fxle->tri=0;
1766 float FXO=0.1f, FXD=0.01f;
1767 QString fxvaluesfile = fxle->resDir.absolutePath()+"_fourxle.fixv";
1768 // qDebug()<<fxvaluesfile ;
1769 if (fxle->keepIso->isChecked()){
1770 if (QFile::exists(fxvaluesfile)){
1771 QFile fxfi(fxvaluesfile);
1772 fxfi.open(QIODevice::ReadOnly|QIODevice::Text);
1773 QStringList lll = QString(fxfi.readAll()).split(' ');
1774 if (lll.size()==2){
1775 FXO=lll.at(0).toFloat();
1776 FXD=lll.at(1).toFloat();
1777 }
1778 fxfi.close();
1779 }else{
1780 QFile fxfi(fxvaluesfile);
1781 fxfi.open(QIODevice::WriteOnly|QIODevice::Text);
1782 FXO = fxle->sigma[1];
1783 FXD = fxle->sigma[0];
1784 fxfi.write(QString("%1 %2").arg(FXO,0,'g').arg(FXD,0,'g').toLocal8Bit());
1785 fxfi.close();
1786 }
1787 }else{
1788 QFile::remove(fxvaluesfile);
1789 }
1790 printf("sigma0 %f sigma1 %f \n",fxle->sigma[0]*2.7f,fxle->sigma[1]*1.2f);
1791 for (int fac=imin; fac<imax; fac++){
1792 switch (fac){
1793 case 0:
1794 if (neu) fxle->iso[1]=(fxle->keepIso->isChecked())?FXD*2.7f:fxle->sigma[0]*2.7f;
1795 else fxle->iso[1]=fabs(fxle->iso[1]);
1796 fxle->mtyp=1;
1797 break;
1798 case 1:
1799 if (neu) fxle->iso[1]=(fxle->keepIso->isChecked())?-FXD*2.7f:-fxle->sigma[0]*2.7f;
1800 else fxle->iso[1]=-fabs(fxle->iso[1]);
1801 fxle->mtyp=1;
1802 break;
1803 case 2:
1804 if (neu) fxle->iso[0]=(fxle->keepIso->isChecked())?FXO*1.2f:fxle->sigma[1]*1.2f;
1805 // printf("blau %g %d\n",fxle->iso[0],neu);
1806 fxle->mtyp=0;
1807 break;
1808 case 3:
1809 if (neu) fxle->iso[0]=(fxle->keepIso->isChecked())?-FXO*1.2f:-fxle->sigma[1]*1.2f;
1810 else fxle->iso[0]=-fabs(fxle->iso[0]);
1811 // printf("orange %g %d\n",fxle->iso[0],neu);
1812 fxle->mtyp=0;
1813 break;
1814 }
1815 {
1816 int ix,iy,iz;
1817 for (iz=0; iz<fxle->n3;iz++){
1818 for (iy=0; iy<fxle->n2;iy++){
1819 for (ix=0; ix<fxle->n1;ix++){
1820 fxle->CalcVertex(ix,iy,iz);
1821 }
1822 }
1823 }
1824 }
1825
1826 int mper=(chgl->niceTrans->isChecked())?6:1;
1827 for (int pers=0;pers<mper;pers++){
1828
1829 if (fxle->nodex==nullptr) {fprintf(stderr,"NODE X is nullptr!\n");return;}
1830 if (fxle->nodey==nullptr) {fprintf(stderr,"NODE Y is nullptr!\n");return;}
1831 if (fxle->nodez==nullptr) {fprintf(stderr,"NODE Z is nullptr!\n");return;}
1832 int facnpers=fac+pers*4;
1833 chgl->fVertexes[facnpers].clear();
1834 chgl->fNormals[facnpers].clear();
1835
1836 int h,k,l;
1837 int step=0;
1838 switch (pers) {
1839 case 0:
1840 //fprintf(stderr, "makeElement case 0 n1=%d n2=%d n3=%d N4=%d fac+pers=%d nodes %p %p %p\n",fxle->n1,fxle->n2,fxle->n3,fxle->n4,facnpers,fxle->nodex,fxle->nodey,fxle->nodez);
1841 for ( h=0; h<fxle->n1;h++){
1842 for ( k=0; k<fxle->n2;k++){
1843 for ( l=0; l<fxle->n3;l++){
1844 if (fxle->nodex==nullptr) {fprintf(stderr,"NODE X is nullptr!\n");return;}
1845 if (fxle->nodey==nullptr) {fprintf(stderr,"NODE Y is nullptr!\n");return;}
1846 if (fxle->nodez==nullptr) {fprintf(stderr,"NODE Z is nullptr!\n");return;}
1847 fxle->MakeElement(h,k,l,facnpers);
1848 }
1849 }
1850 step++;step%=mper;
1851 }
1852 //printf("makeElement made %d Vertexes %p %p %p\n", chgl->fVertexes[facnpers].size(), fxle->nodex, fxle->nodey, fxle->nodez);
1853 break;
1854 case 1:
1855 //fprintf(stderr, "makeElement case 1\n");
1856 for ( h=fxle->n1; h>=0;h--){
1857 for ( k=fxle->n2; k>=0;k--){
1858 for ( l=fxle->n3; l>=0;l--){
1859 if (fxle->nodex==nullptr) {fprintf(stderr,"NODE X is nullptr!\n");return;}
1860 if (fxle->nodey==nullptr) {fprintf(stderr,"NODE Y is nullptr!\n");return;}
1861 if (fxle->nodez==nullptr) {fprintf(stderr,"NODE Z is nullptr!\n");return;}
1862 fxle->MakeElement(h,k,l,facnpers);
1863 }
1864 }
1865 step++;step%=mper;
1866 }
1867 break;
1868 case 2:
1869 //fprintf(stderr, "makeElement case 2\n");
1870 for ( k=0; k<fxle->n2;k++){
1871 for ( h=0; h<fxle->n1;h++){
1872 for ( l=0; l<fxle->n3;l++){
1873 if (fxle->nodex==nullptr) {fprintf(stderr,"NODE X is nullptr!\n");return;}
1874 if (fxle->nodey==nullptr) {fprintf(stderr,"NODE Y is nullptr!\n");return;}
1875 if (fxle->nodez==nullptr) {fprintf(stderr,"NODE Z is nullptr!\n");return;}
1876 fxle->MakeElement(h,k,l,facnpers);
1877 }
1878 }
1879 step++;step%=mper;
1880 }
1881 break;
1882 case 3:
1883 //fprintf(stderr, "makeElement case 3\n");
1884 for ( k=fxle->n2; k>=0;k--){
1885 for ( h=fxle->n1; h>=0;h--){
1886 for ( l=fxle->n3; l>=0;l--){
1887 if (fxle->nodex==nullptr) {fprintf(stderr,"NODE X is nullptr!\n");return;}
1888 if (fxle->nodey==nullptr) {fprintf(stderr,"NODE Y is nullptr!\n");return;}
1889 if (fxle->nodez==nullptr) {fprintf(stderr,"NODE Z is nullptr!\n");return;}
1890 fxle->MakeElement(h,k,l,facnpers);
1891 }
1892 }
1893 step++;step%=mper;
1894 }
1895 break;
1896 case 4:
1897 //fprintf(stderr, "makeElement case 4\n");
1898 for ( l=0; l<fxle->n3;l++){
1899 for ( h=0; h<fxle->n1;h++){
1900 for ( k=0; k<fxle->n2;k++){
1901 if (fxle->nodex==nullptr) {fprintf(stderr,"NODE X is nullptr!\n");return;}
1902 if (fxle->nodey==nullptr) {fprintf(stderr,"NODE Y is nullptr!\n");return;}
1903 if (fxle->nodez==nullptr) {fprintf(stderr,"NODE Z is nullptr!\n");return;}
1904 fxle->MakeElement(h,k,l,facnpers);
1905 }
1906 }
1907 step++;step%=mper;
1908 }
1909 break;
1910 case 5:
1911 //fprintf(stderr, "makeElement case 5\n");
1912 for ( l=fxle->n3; l>=0;l--){
1913 for ( h=fxle->n1; h>=0;h--){
1914 for ( k=fxle->n2; k>=0;k--){
1915 if (fxle->nodex==nullptr) {fprintf(stderr,"NODE X is nullptr!\n");return;}
1916 if (fxle->nodey==nullptr) {fprintf(stderr,"NODE Y is nullptr!\n");return;}
1917 if (fxle->nodez==nullptr) {fprintf(stderr,"NODE Z is nullptr!\n");return;}
1918 fxle->MakeElement(h,k,l,facnpers);
1919 }
1920 }
1921 step++;step%=mper;
1922 }
1923 break;
1924 }
1925 }
1926 }
1927 fxle->iso[0]=fabs(fxle->iso[0]);
1928
1929 fxle->info=QString("<b>%13 Map:</b><font color=\"%12\">%1 eÅ<sup>-3</sup></font>"
1930 "<font color=\"%11\"> %2 eÅ<sup>-3</sup> (σ = %7) </font><br><font color=grey> Hint: [%3 Scroll (up or down)] to change. </font><br>"
1931 "<b>Fo-Map:</b><font color=\"%10\">%4 eÅ<sup>-3</sup> (σ = %8)</font><br><font color=grey> Hint: [%5 Scroll (up or down)] to change. </font><br>Time for creating map surfaces <b>%6 s</b>. %9 Triangles drawn.<br>")
1932 .arg(fxle->iso[1],6,'g',2)
1933 .arg(-fxle->iso[1],6,'g',2)
1934 .arg(QKeySequence(Qt::ControlModifier).toString(QKeySequence::NativeText))
1935 .arg(fxle->iso[0],6,'g',2)
1936 .arg(QKeySequence(Qt::ShiftModifier).toString(QKeySequence::NativeText))
1937 .arg(fxle->suti.elapsed()/1000.0,1,'f',1)
1938 .arg(fxle->sigma[0],6,'g',2)
1939 .arg(fxle->sigma[1],6,'g',2)
1940 .arg(fxle->tri)
1941 .arg(chgl->fopc.name())
1942 .arg(chgl->dipc.name())
1943 .arg(chgl->dimc.name())
1944 .arg((chgl->dipc.name() == QString("#008080"))?"BODD":"Fo-Fc1")//quick and dirty check for BODD
1945 ;
1946 //chgl->pause=false;
1947 fxle->connect(chgl,SIGNAL(inimibas()),fxle,SLOT(inimap()));
1948 fxle->connect(chgl,SIGNAL(neuemitte(V3)),fxle, SLOT(bewegt(V3)));
1949 fxle->connect(chgl,SIGNAL(diffscroll(int ,int )),fxle,SLOT(change_iso(int ,int )));
1950 chgl->fofcact->setVisible(!chgl->fVertexes[0].isEmpty()||!chgl->fVertexes[1].isEmpty());
1951 chgl->foact->setVisible(!chgl->fVertexes[2].isEmpty());
1952 //printf("wer ist empty %d %d %d %d\n",!chgl->fVertexes[0].isEmpty(),!chgl->fVertexes[1].isEmpty(),!chgl->fVertexes[2].isEmpty(),!chgl->fVertexes[0].isEmpty()||!chgl->fVertexes[1].isEmpty());
1953 // printf("Fomap ist da? %d Fo-Fcmap ist da? %d \n",chgl->fVertexes[2].isEmpty(),chgl->fVertexes[0].isEmpty());
1954
1955 fxle->acnt=chgl->objCnt;
1956 }
1957
1958
reportInfo()1959 void FourXle::reportInfo(){
1960 emit bigmessage(info);
1961 }
1962
change_iso(int numsteps,int diff)1963 void FourXle::change_iso(int numsteps,int diff){
1964 /*! canges the iso value of the fo or fo-fc maps and redraws them
1965 */
1966 //printf("isochaged\n");
1967 iso[diff]=fabs(iso[diff]);
1968 iso[diff]+=iso[diff]*numsteps/10.0f;
1969 int mi=0,ma=5;
1970 switch (diff) {
1971 case 0:
1972 mi=2;ma=4;
1973 break;
1974 case 1:
1975 mi=0; ma=2;
1976 break;
1977
1978 }
1979 gen_surface(false,mi,ma);
1980 chgl->updateGL();
1981 }
1982
foactDestroyed()1983 void FourXle::foactDestroyed() {
1984 if (chgl!=nullptr) chgl->foact = nullptr;
1985 }
1986
fofcactDestroyed()1987 void FourXle::fofcactDestroyed() {
1988 if (chgl!=nullptr)chgl->fofcact = nullptr;
1989 }
1990
chglDestroyed()1991 void FourXle::chglDestroyed() {
1992 //fprintf(stderr,"chgl is dead omg ...argh...\n");
1993
1994 chgl = nullptr;
1995 }
1996
CalcVertex(int ix,int iy,int iz)1997 void FourXle::CalcVertex( int ix, int iy, int iz) {
1998 V3 mdz=(0.5*dx)+(0.5*dy)+(0.5*dz);
1999 V3 o,fl,m2u=V3(0.5,0.5,0.5);
2000 mole->frac2kart(m2u,m2u);
2001 double vo=0, vx=0,vy=0,vz=0;
2002 int idx=dex(ix,iy,iz);
2003 nodex[idx].flag=0;
2004 nodey[idx].flag=0;
2005 nodez[idx].flag=0;
2006 if (mtyp==0){//*datfo,*datfo_fc,*datf1_f2
2007 vo = datfo[ idx] -iso[mtyp];
2008 vx = datfo[ dex(ix+1,iy,iz)] -iso[mtyp];
2009 vy = datfo[ dex(ix,iy+1,iz)] -iso[mtyp];
2010 vz = datfo[ dex(ix,iy,iz+1)] -iso[mtyp];
2011 }else {
2012 vo = datfo_fc[idx] -iso[mtyp];
2013 vx = datfo_fc[dex(ix+1,iy,iz)] -iso[mtyp];
2014 vy = datfo_fc[dex(ix,iy+1,iz)] -iso[mtyp];
2015 vz = datfo_fc[dex(ix,iy,iz+1)] -iso[mtyp];
2016 }
2017 V3 nor=V3(0,0,0);//Normalize((vx-vo)*dx+(vy-vo)*dy+(vz-vo)*dz);
2018 if (Intersect(vo,vx)) {
2019 if (maptrunc==2) o=dx*((vo/(vo-vx))+ix) + dy*iy + dz*iz+m2u;
2020 else o=dx*((vo/(vo-vx))+ix) + dy*iy + dz*iz+urs;
2021
2022 mole->kart2frac(o,o);
2023 o+=V3(-0.5,-0.5,-0.5);
2024 fl=V3(floor(o.x),floor(o.y),floor(o.z));
2025 o+=-1.0*fl;
2026 o+=V3(0.5,0.5,0.5);
2027 mole->frac2kart(o,o);
2028 if (maptrunc!=2) o+=-1.0*urs;
2029 else o+=-1.0*m2u;
2030 o+=mdz;
2031 // orte.append(o);
2032 nodex[idx].vertex=o;
2033 nodex[idx].normal=nor;
2034 nodex[idx].flag=1;
2035 }
2036 if (Intersect(vo,vy)) {
2037 if (maptrunc==2)o=dx*ix + dy*((vo/(vo-vy))+iy) + dz*iz+m2u;
2038 else o=dx*ix + dy*((vo/(vo-vy))+iy) + dz*iz+urs;
2039 mole->kart2frac(o,o);
2040 o+=V3(-0.5,-0.5,-0.5);
2041 fl=V3(floor(o.x),floor(o.y),floor(o.z));
2042 o+=-1.0*fl;
2043 o+=V3(0.5,0.5,0.5);
2044 mole->frac2kart(o,o);
2045 if (maptrunc!=2) o+=-1.0*urs;
2046 else o+=-1.0*m2u;
2047 o+=mdz;
2048 // orte.append(o);
2049 nodey[idx].vertex=o;
2050 nodey[idx].normal=nor;
2051 nodey[idx].flag=1;
2052 }
2053 if (Intersect(vo,vz)) {
2054 if (maptrunc==2)o=dx*ix + dy*iy + dz*((vo/(vo-vz))+iz)+m2u;
2055 else o=dx*ix + dy*iy + dz*((vo/(vo-vz))+iz)+urs;
2056 mole->kart2frac(o,o);
2057 o+=V3(-0.5,-0.5,-0.5);
2058 fl=V3(floor(o.x),floor(o.y),floor(o.z));
2059 o+=-1.0*fl;
2060 o+=V3(0.5,0.5,0.5);
2061 mole->frac2kart(o,o);
2062 if (maptrunc!=2) o+=-1.0*urs;
2063 else o+=-1.0*m2u;
2064 o+=mdz;
2065 // orte.append(o);
2066 nodez[idx].vertex=o;
2067 nodez[idx].normal=nor;
2068 nodez[idx].flag=1;
2069 }
2070 }
2071
makeFaces(int n,int faps,FNode poly[])2072 void FourXle::makeFaces(int n, int faps, FNode poly[] ){
2073 if (n<3) return; //less then 3 vertices -> nothing to do
2074 V3 mid_ver = V3(0,0,0);
2075 V3 mid_nor2 = V3(0,0,0);
2076 for(int i=0; i<n; i++ ){
2077 if ((i>0)&&(i<n-1)) mid_nor2 += Normalize((poly[i].vertex-poly[0].vertex)%(poly[i+1].vertex-poly[0].vertex));
2078 mid_ver += poly[i].vertex;
2079 }
2080 mid_nor2=Normalize(mid_nor2);
2081 for(int i=0; i<n; i++ ){
2082 poly[i].normal = mid_nor2;
2083 }
2084 mid_ver *= (1.0/n);
2085 //mid_nor2= Normalize(mid_nor2);
2086 V3 mit=V3(1,1,1);
2087 mole->frac2kart(mit,mit);
2088 for (int w=0; w<27; w++){
2089 if (maptrunc==0){w=13;}
2090 else if ((maptrunc==1)&&(Distance(mit-urs,mid_ver+delDA[w])>(map_radius*map_radius))) continue;
2091 else if (maptrunc==2) {
2092 mole->kart2frac(mid_ver+delDA[w],oc);
2093 int ax = static_cast<int>((oc.x)/dxc.x-0.499),
2094 bx = static_cast<int>((oc.y)/dyc.y-0.499),
2095 cx = static_cast<int>((oc.z)/dzc.z-0.499);
2096 int DEX=dex3(ax,bx,cx);
2097 if (!noGoMap[DEX]) continue;
2098 }
2099 for (int k=1; k<n-1;k++){
2100 tri++;
2101
2102 chgl->fVertexes[faps].append(poly[k].vertex.x+delDA[w].x);
2103 chgl->fVertexes[faps].append(poly[k].vertex.y+delDA[w].y);
2104 chgl->fVertexes[faps].append(poly[k].vertex.z+delDA[w].z);
2105
2106 chgl->fVertexes[faps].append(poly[(k+1)].vertex.x+delDA[w].x);
2107 chgl->fVertexes[faps].append(poly[(k+1)].vertex.y+delDA[w].y);
2108 chgl->fVertexes[faps].append(poly[(k+1)].vertex.z+delDA[w].z);
2109
2110 chgl->fVertexes[faps].append(poly[0].vertex.x+delDA[w].x);
2111 chgl->fVertexes[faps].append(poly[0].vertex.y+delDA[w].y);
2112 chgl->fVertexes[faps].append(poly[0].vertex.z+delDA[w].z);
2113
2114
2115 chgl->fNormals[faps].append(poly[k].normal.x);
2116 chgl->fNormals[faps].append(poly[k].normal.y);
2117 chgl->fNormals[faps].append(poly[k].normal.z);
2118
2119 chgl->fNormals[faps].append(poly[(k+1)].normal.x);
2120 chgl->fNormals[faps].append(poly[(k+1)].normal.y);
2121 chgl->fNormals[faps].append(poly[(k+1)].normal.z);
2122
2123 chgl->fNormals[faps].append(poly[0].normal.x);
2124 chgl->fNormals[faps].append(poly[0].normal.y);
2125 chgl->fNormals[faps].append(poly[0].normal.z);
2126 }
2127 if (maptrunc==0) w=27;
2128 }
2129 }
2130
IndexSelected(FNode & node0,FNode & node1,FNode & node2,FNode & node3)2131 int FourXle::IndexSelected( FNode& node0, FNode& node1, FNode& node2, FNode& node3 ) {
2132 if( node1 && node2 && node3 ){
2133 double d1 = Distance( node0.vertex, node1.vertex) + Distance( node3.vertex, node2.vertex) ;
2134 double d2 = Distance( node0.vertex, node2.vertex) + Distance( node3.vertex, node1.vertex );
2135 if( d1 > d2 ) return 2; else return 1;
2136 }else{
2137 if( node1 ) return 1;
2138 else if( node2 ) return 2;
2139 else if( node3 ) return 3;
2140 }
2141 return 0;
2142
2143 }
2144
MakeElement(int ix,int iy,int iz,int faps)2145 void FourXle::MakeElement(int ix, int iy, int iz , int faps) {
2146 QMutexLocker locker(&mutex);
2147 int conn[12][2][4] = {
2148 {{ 0, 1, 7, 6}, { 0, 2, 8, 3}}, // 0
2149 {{ 1, 2, 5, 4}, { 1, 0, 6, 7}}, // 1
2150 {{ 2, 0, 3, 8}, { 2, 1, 4, 5}}, // 2
2151 {{ 3, 8, 2, 0}, { 3, 4,10, 9}}, // 3
2152 {{ 4, 3, 9,10}, { 4, 5, 2, 1}}, // 4
2153 {{ 5, 4, 1, 2}, { 5, 6, 9,11}}, // 5
2154 {{ 6, 5,11, 9}, { 6, 7, 1, 0}}, // 6
2155 {{ 7, 6, 0, 1}, { 7, 8,11,10}}, // 7
2156 {{ 8, 7,10,11}, { 8, 3, 0, 2}}, // 8
2157 {{ 9,10, 4, 3}, { 9,11, 5, 6}}, // 9
2158 {{10,11, 8, 7}, {10, 9, 3, 4}}, // 10
2159 {{11, 9, 6, 5}, {11,10, 7, 8}} // 11
2160 };
2161 FNode node[12];
2162 FNode polygon[12];
2163 /*
2164 if (nodex==nullptr) {fprintf(stderr,"NODE X is nullptr!\n");return;}
2165 if (nodey==nullptr) {fprintf(stderr,"NODE Y is nullptr!\n");return;}
2166 if (nodez==nullptr) {fprintf(stderr,"NODE Z is nullptr!\n");return;}
2167 if (
2168 ((ix+iy*s1+iz*n4) %n5)<0||
2169 ((ix+iy*s1+iz*n4) %n5)<0||
2170 ((ix+iy*s1+iz*n4) %n5)<0||
2171 ((ix+iy*s1+((iz+1)%n3)*n4) %n5)<0||
2172 ((ix+iy*s1+((iz+1)%n3)*n4) %n5)<0||
2173 ((ix+((iy+1)%n2)*s1+iz*n4) %n5)<0||
2174 ((ix+((iy+1)%n2)*s1+iz*n4) %n5)<0||
2175 ((((1+ix)%n1)+iy*s1+iz*n4) %n5)<0||
2176 ((((1+ix)%n1)+iy*s1+iz*n4) %n5)<0||
2177 ((ix+((iy+1)%n2)*s1+((iz+1)%n3)*n4)%n5)<0||
2178 ((((ix+1)%n1)+iy*s1+((iz+1)%n3)*n4)%n5)<0||
2179 ((((ix+1)%n1)+((iy+1)%n2)*s1+iz*n4)%n5)<0
2180 )printf("EEEEK!\n");
2181 */
2182 node[ 0] = nodex[(ix+iy*n1+iz*n4) %n5]; // 000x
2183 node[ 1] = nodey[(ix+iy*n1+iz*n4) %n5]; // 000y
2184 node[ 2] = nodez[(ix+iy*n1+iz*n4) %n5]; // 000z
2185 node[ 3] = nodex[(ix+iy*n1+((iz+1)%n3)*n4) %n5]; // 001y
2186 node[ 4] = nodey[(ix+iy*n1+((iz+1)%n3)*n4) %n5]; // 001z
2187 node[ 5] = nodez[(ix+((iy+1)%n2)*n1+iz*n4) %n5]; // 010x
2188 node[ 6] = nodex[(ix+((iy+1)%n2)*n1+iz*n4) %n5]; // 010y
2189 node[ 7] = nodey[(((1+ix)%n1)+iy*n1+iz*n4) %n5]; // 100y
2190 node[ 8] = nodez[(((1+ix)%n1)+iy*n1+iz*n4) %n5]; // 100z
2191 node[ 9] = nodex[(ix+((iy+1)%n2)*n1+((iz+1)%n3)*n4)%n5];// 011x
2192 node[10] = nodey[(((ix+1)%n1)+iy*n1+((iz+1)%n3)*n4)%n5];// 101y
2193 node[11] = nodez[(((ix+1)%n1)+((iy+1)%n2)*n1+iz*n4)%n5];// 110z
2194 if ((static_cast<char>(node[0])+node[1]+node[2]+node[3]+node[4]+node[5]
2195 +node[6]+node[7]+node[8]+node[9]+node[10]+node[11])==0) return;
2196 for( int is=0; is<12; is++ ) {
2197 if( !node[is] ) continue;
2198
2199 int n=0, i=is, m=0;//,ai=i;
2200 double dis;
2201 dis=0.0;
2202 do {
2203 polygon[n++]= node[i];
2204 int sol = IndexSelected(
2205 node[conn[i][m][0]],
2206 node[conn[i][m][1]],
2207 node[conn[i][m][2]],
2208 node[conn[i][m][3]]);
2209 //ai=i;
2210 i = conn[i][m][sol];
2211 if( sol == 2 ) m ^= 1;
2212 dis+=Distance(polygon[0].vertex,node[i].vertex);
2213 node[i].flag= 0;
2214 }while( (i!=is)&&(n<11) );
2215 if (n>=3) {
2216 if (dis<5){
2217 //for (int polni=0; polni<n; polni++){polygon[polni].normal+=Normalize((polygon[(polni+1)%n].vertex-polygon[polni].vertex)%(polygon[(polni+n-1)%n].vertex-polygon[polni].vertex));}
2218 makeFaces( n, faps, polygon );
2219 }//*
2220 else {
2221 int axe=0;
2222 double delx=0,dely=0,delz=0;
2223 double mind=100000000;
2224 V3 minp=V3(10000,10000,10000),lihiun=V3(-10000,-10000,-10000);
2225 int minii=0;
2226 for (int polni=1; polni<=n; polni++){
2227 delx+=fabs(polygon[polni-1].vertex.x - polygon[polni%n].vertex.x);
2228 dely+=fabs(polygon[polni-1].vertex.y - polygon[polni%n].vertex.y);
2229 delz+=fabs(polygon[polni-1].vertex.z - polygon[polni%n].vertex.z);
2230 if (Distance(polygon[polni%n].vertex,lihiun)<mind) {
2231 mind=Distance(polygon[polni%n].vertex,minp);
2232 minii=polni%n;
2233 }
2234 }
2235 minp=polygon[minii].vertex;
2236 axe|=(delx>1)?1:0;
2237 axe|=(dely>1)?2:0;
2238 axe|=(delz>1)?4:0;
2239 for (int polni=0; polni<n; polni++){
2240 V3 neo=polygon[polni].vertex;
2241 double lang =Distance(minp,neo);
2242 if ((lang>Distance(minp,neo+dx*n1))&&(axe&1)) neo+=dx*n1;
2243 else if ((lang>Distance(minp,neo-dx*n1))&&(axe&1)) neo+=-n1*dx;
2244 lang =Distance(minp,neo);
2245 if ((lang>Distance(minp,neo+dy*n2))&&(axe&2)) neo+=n2*dy;
2246 else if ((lang>Distance(minp,neo-dy*n2))&&(axe&2)) neo+=-n2*dy;
2247 lang =Distance(minp,neo);
2248 if ((lang>Distance(minp,neo+n3*dz))&&(axe&4)) neo+=n3*dz;
2249 else if ((lang>Distance(minp,neo-n3*dz))&&(axe&4)) neo+=-n3*dz;
2250 polygon[polni].vertex=neo;
2251 }
2252 dis=0;
2253 for (int polni=1; polni<=n; polni++){
2254 dis+=Distance(polygon[polni-1].vertex , polygon[polni%n].vertex);
2255 }
2256 if (dis<5)
2257 makeFaces( n, faps, polygon );
2258 }
2259 // */
2260 }
2261 }
2262
2263 }
2264
sigmaFo()2265 float FourXle::sigmaFo(){
2266 float DM=0.0f;
2267 float DS=0.0f;
2268 float f=0.0f;
2269 for (int i=0; i<n5; i++){
2270 f=datfo[i];
2271 DM+=f;
2272 DS+=f*f;
2273 }
2274 return sqrt((DS/n5)-((DM/n5)*(DM/n5)));
2275 }
2276
sigmaFoFc()2277 float FourXle::sigmaFoFc(){
2278 float DM=0.0f;
2279 float DS=0.0f;
2280 float f=0.0f;
2281 for (int i=0; i<n5; i++){
2282 f=datfo_fc[i];
2283 DM+=f;
2284 DS+=f*f;
2285 }
2286 return sqrt((DS/n5)-((DM/n5)*(DM/n5)));
2287 }
2288
jnk()2289 void FourXle::jnk(){
2290 if(!n5)return;
2291 float mini=9e37f,maxi=-9e37f;
2292 float f,r,fstep;
2293 double DM,DS,sigma,w,e_net,e_gross;
2294 r=powf((3*(n1-1)*(n2-1)*(n3-1)),1.0f/3.0f);
2295 // printf("%d %d %d %d %f\n",n1,n2,n3,n5,r);
2296 float *datfo_fcstp=(float*) malloc(sizeof(float)*n5);
2297 if (datfo_fcstp==nullptr)return ;
2298 float invstep=100.0f;
2299 float step=1.0f/invstep;
2300 float rhomind2=99999.0f,rhomaxd2=-99999.0f;
2301 QMap<float,int> hash;
2302 QMap<float,float> hashf;
2303 DM=0.0;
2304 DS=0.0;
2305 w=0.0;
2306 for (int i=0; i<n5; i++){
2307 f=datfo_fc[i];
2308 mini=qMin(f,mini);
2309 maxi=qMax(f,maxi);
2310 DM+=f;
2311 DS+=f*f;
2312 w+=fabs(f);
2313 datfo_fcstp[i]=floorf(f*invstep);
2314 }// */
2315 e_net=(DM/n5)*C[14];//C[14] ist das Volumen der UC
2316 e_gross=w/(2*n5)*C[14];
2317 sigma=sqrt((DS/n5)-((DM/n5)*(DM/n5)));
2318 //printf("%g %g %g",e_net,e_gross,sigma);
2319 for (int zi=0;zi<n3;zi++)
2320 for (int yi=0;yi<n2;yi++)
2321 for (int xi=0;xi<(n1-1);xi++){
2322 if (datfo_fcstp[dex(xi,yi,zi)]>datfo_fcstp[dex(xi+1,yi,zi)]){
2323 fstep=datfo_fcstp[dex(xi,yi,zi)];
2324 int ze=fstep-datfo_fcstp[dex(xi+1,yi,zi)];
2325 for (int zii=0; zii<ze; zii++) hash[(fstep-zii)*step]++;
2326 }else if (datfo_fcstp[dex(xi,yi,zi)]<datfo_fcstp[dex(xi+1,yi,zi)]){
2327 fstep=datfo_fcstp[dex(xi+1,yi,zi)];
2328 int ze=fstep-datfo_fcstp[dex(xi,yi,zi)];
2329 for (int zii=0; zii<ze; zii++) hash[(fstep-zii)*step]++;
2330 }
2331 }
2332 for (int xi=0;xi<n1;xi++)
2333 for (int zi=0;zi<n3;zi++)
2334 for (int yi=0;yi<(n2-1);yi++){
2335 if (datfo_fcstp[dex(xi,yi,zi)]>datfo_fcstp[dex(xi,yi+1,zi)]){
2336 fstep=datfo_fcstp[dex(xi,yi,zi)];
2337 int ze=fstep-datfo_fcstp[dex(xi,yi+1,zi)];
2338 for (int zii=0; zii<ze; zii++) hash[(fstep-zii)*step]++;
2339 }else if (datfo_fcstp[dex(xi,yi,zi)]<datfo_fcstp[dex(xi,yi+1,zi)]){
2340 fstep=datfo_fcstp[dex(xi,yi+1,zi)];
2341 int ze=fstep-datfo_fcstp[dex(xi,yi,zi)];
2342 for (int zii=0; zii<ze; zii++) hash[(fstep-zii)*step]++;
2343 }
2344 }
2345
2346 for (int yi=0;yi<n2;yi++)
2347 for (int xi=0;xi<n1;xi++)
2348 for (int zi=0;zi<(n3-1);zi++){
2349 if (datfo_fcstp[dex(xi,yi,zi)]>datfo_fcstp[dex(xi,yi,zi+1)]){
2350 fstep=datfo_fcstp[dex(xi,yi,zi)];
2351 int ze=fstep-datfo_fcstp[dex(xi,yi,zi+1)];
2352 for (int zii=0; zii<ze; zii++) hash[(fstep-zii)*step]++;
2353 }else if (datfo_fcstp[dex(xi,yi,zi)]<datfo_fcstp[dex(xi,yi,zi+1)]){
2354 fstep=datfo_fcstp[dex(xi,yi,zi+1)];
2355 int ze=fstep-datfo_fcstp[dex(xi,yi,zi)];
2356 for (int zii=0; zii<ze; zii++) hash[(fstep-zii)*step]++;
2357 }
2358 }
2359 free(datfo_fcstp);
2360 //float range=fmaxf(fabsf(mini),fabsf(maxi));
2361 QDialog *jnkdlg = new QDialog();
2362 QGraphicsScene*scene= new QGraphicsScene(-30,-50,550,586);
2363 scene->setBackgroundBrush(QBrush(QColor("#e9f7d6")));
2364 scene->clear ();
2365 QGraphicsItem *itm;
2366 for (int i=0; i<21;i++){
2367 itm=scene->addLine(i*25,0,i*25,500,(i%5)?QPen(QColor("#cbdbbb"),0):QPen(QColor("#959d9d"),0));
2368 itm->setData(0,-1);
2369 }
2370 for (int i=0; i<31;i++){
2371 itm=scene->addLine(0,i*16.66666666666667,500,i*16.66666666666667,(i%5)?QPen(QColor("#cbdbbb"),0):QPen(QColor("#959d9d"),0));
2372 itm->setData(0,-1);
2373 }
2374 itm=scene->addLine(250,0,250,500,QPen(QColor("#000000"),0));
2375 itm=scene->addLine(0,83.33333333333333,500,83.33333333333333,QPen(QColor("#000000"),0));
2376 itm=scene->addLine(0,250,500,250,QPen(QColor("#000000"),0));
2377 itm=scene->addLine(0,416.6666666666667,500,416.6666666666667,QPen(QColor("#000000"),0));
2378 QGraphicsTextItem *txt = scene->addText("3");
2379 txt->moveBy(-15,73.33333333333333);
2380 txt = scene->addText("2");
2381 txt->moveBy(-15,240);
2382 txt = scene->addText("1");
2383 txt->moveBy(-15,406.6666666666667);
2384 txt = scene->addText("-1.0");
2385 txt->moveBy(-15,500);
2386 txt = scene->addText("0");
2387 txt->moveBy(241,500);
2388 txt = scene->addText("1.0");
2389 txt->moveBy(485,500);
2390
2391
2392 QMapIterator<float, int> i(hash);
2393 while (i.hasNext()) {
2394 i.next();
2395 f=hashf[i.key()]=logf(i.value())/logf(r);
2396 scene->addEllipse(250+(i.key())*250,500-((f-0.5)/3.)*500,4,4,QPen(Qt::NoPen),QBrush(QColor("#0907e6")));
2397
2398 if (f>2){
2399 rhomind2=qMin(rhomind2,i.key());
2400 rhomaxd2=qMax(rhomaxd2,i.key());
2401 }
2402 // printf("jnk: %g %g %g %g %g %g\n",i.key(),f,rhomind2,rhomaxd2,250+(i.key())*250,500-((f-0.5)/3.)*500);
2403 }
2404 // printf("\ndf(0)= %g %d %g %g\n",hashf.value(0.0f),hash.value(0.0f),rhomind2,rhomaxd2);
2405
2406 float m,b;
2407 m=(hashf.value(rhomind2)-hashf.value(rhomind2-step))/step;
2408 b=hashf.value(rhomind2)-m*rhomind2;
2409 rhomind2=(2.0f-b)/m;
2410 m=(hashf.value(rhomaxd2)-hashf.value(rhomaxd2-step))/step;
2411 b=hashf.value(rhomaxd2)-m*(rhomaxd2);
2412 rhomaxd2=(2.0f-b)/m;
2413 // printf("df(0)= %g %d %g %g\n",hashf.value(0.0f),hash.value(0.0f),rhomind2,rhomaxd2);
2414 txt = scene->addText(QString("df(0) = %1").arg((double)hashf.value(0.0f),7,'f',2),QFont("Courier",10));
2415 txt->moveBy(328,2);
2416 txt = scene->addText(QString("min(d=2) = %1 eA^-3").arg((double)rhomind2,7,'f',3),QFont("Courier",10));
2417 txt->moveBy(305,18);
2418 txt = scene->addText(QString("max(d=2) = %1 eA^-3").arg((double)rhomaxd2,7,'f',3),QFont("Courier",10));
2419 txt->moveBy(305,34);
2420 txt = scene->addText(QString("min = %1 eA^-3").arg((double)mini,7,'f',3),QFont("Courier",10));
2421 txt->moveBy(305,50);
2422 txt = scene->addText(QString("max = %1 eA^-3").arg((double)maxi,7,'f',3),QFont("Courier",10));
2423 txt->moveBy(305,66);
2424 txt = scene->addText(QString("e_net = %1 e").arg((double)e_net,8,'g',3),QFont("Courier",10));
2425 txt->moveBy(295,82);
2426 txt = scene->addText(QString("e_gross = %1 e").arg((double)e_gross,8,'g',3),QFont("Courier",10));
2427 txt->moveBy(295,98);
2428 txt = scene->addText(QString(" = %1 eA^-3").arg((double)sigma,8,'f',3),QFont("Courier",10));
2429 txt->moveBy(295,116);
2430 txt = scene->addText(QString("nx ny nz %1 x %2 x %3 = %4").arg(n1).arg(n2).arg(n3).arg(n5),QFont("Courier",9));
2431 txt->moveBy(30,2);
2432 txt = scene->addText(QString("fractal dimension vs. residual density"),QFont("Helvetica",14,QFont::Bold));
2433 txt->moveBy(100,-40);
2434 txt = scene->addText(QString("df"),QFont("Helvetica",12,QFont::Bold));
2435 txt->moveBy(-20,20);
2436 txt = scene->addText(QString("0"),QFont("Helvetica",12,QFont::Bold));
2437 txt->moveBy(460,500);
2438
2439 #ifdef __APPLE__
2440 txt = scene->addText(QString::fromUtf8("��"),QFont("Helvetica",16,QFont::Bold));
2441 txt->moveBy(450,490);
2442 txt = scene->addText(QString::fromUtf8("��"),QFont("Helvetica",14));
2443 txt->moveBy(295,8);
2444 txt = scene->addText(QString::fromUtf8("��"),QFont("Helvetica",14));
2445 txt->moveBy(295,24);
2446 txt = scene->addText(QString::fromUtf8("��"),QFont("Helvetica",14));
2447 txt->moveBy(295,40);
2448 txt = scene->addText(QString::fromUtf8("��"),QFont("Helvetica",14));
2449 txt->moveBy(295,56);
2450 txt = scene->addText(QString::fromUtf8("��"),QFont("Helvetica",14));
2451 txt->moveBy(295,106);
2452 #else
2453 txt = scene->addText(QString::fromUtf8("ρ"),QFont("Helvetica",16,QFont::Bold));
2454 txt->moveBy(450,490);
2455 txt = scene->addText(QString::fromUtf8("ρ"),QFont("Helvetica",14));
2456 txt->moveBy(295,8);
2457 txt = scene->addText(QString::fromUtf8("ρ"),QFont("Helvetica",14));
2458 txt->moveBy(295,24);
2459 txt = scene->addText(QString::fromUtf8("ρ"),QFont("Helvetica",14));
2460 txt->moveBy(295,40);
2461 txt = scene->addText(QString::fromUtf8("ρ"),QFont("Helvetica",14));
2462 txt->moveBy(295,56);
2463 txt = scene->addText(QString::fromUtf8("σ"),QFont("Helvetica",14));
2464 txt->moveBy(295,106);
2465 #endif // */
2466 txt = scene->addText(QString("Please cite as: 'K. Meindl, J. Henn, Acta Cryst., 2008, A64, 404-418.'"),QFont("Helvetica",9));
2467 txt->moveBy(100,516);
2468 QGraphicsView *view = new QGraphicsView(scene,jnkdlg);
2469
2470 QVBoxLayout *lt = new QVBoxLayout(jnkdlg);
2471 lt->addWidget(view);
2472 jnkdlg->show();
2473
2474
2475 }
2476
2477
dex(int x,int y,int z)2478 inline int FourXle::dex(int x,int y, int z){
2479 /*! dex is used to adress elemennts of a one dimensional array by three indizes like it is a 3 dimensional array
2480 * @param x,y,z tree dimensional indices
2481 * if x is < 0 or > n1 it is not a problem because % is used to clamp it.
2482 * if y is < 0 or > n2 it is not a problem because % is used to clamp it.
2483 * if z is < 0 or > n3 it is not a problem because % is used to clamp it.
2484 * \returns index of an 1 dimensional array
2485 */
2486 x=(x+n1)%n1;
2487 y=(y+n2)%n2;
2488 z=(z+n3)%n3;
2489 return x+n1*(y+n2*z);
2490 }
2491
dex3(int x,int y,int z)2492 inline int FourXle::dex3(int x,int y, int z){
2493 int n31=3*n1,n32=3*n2,n33=3*n3;
2494 x=(x+n31)%n31;
2495 y=(y+n32)%n32;
2496 z=(z+n33)%n33;
2497 return x+n31*(y+n32*z);
2498 }
2499
2500 /*=========================P=L=A=N=E===========================================/
2501 *
2502 * Contours lines on a plane
2503 *
2504 *=============================================================================*/
2505
IndexSelectedP(Node & node0,Node & node1,Node & node2,Node & node3)2506 int FourXle::IndexSelectedP( Node& node0, Node& node1, Node& node2, Node& node3 ) {
2507 if( node1 && node2 && node3 ){
2508 GLfloat d1 = Distance( Planorte.at(node0.index).vertex, Planorte.at(node1.index).vertex ) +
2509 Distance( Planorte.at(node3.index).vertex, Planorte.at(node2.index).vertex );
2510 GLfloat d2 = Distance( Planorte.at(node0.index).vertex, Planorte.at(node2.index).vertex ) +
2511 Distance( Planorte.at(node3.index).vertex, Planorte.at(node1.index).vertex );
2512 if( d1 > d2 ) return 2; else return 1;
2513 }else{
2514 if( node1 ) return 1;
2515 else if( node2 ) return 2;
2516 else if( node3 ) return 3;
2517 }
2518 return 0;
2519 }
2520
cutTriangle(GLfloat va,GLfloat vb,GLfloat vc)2521 int cutTriangle(GLfloat va, GLfloat vb, GLfloat vc){
2522 if ((va<=0.0f)&&(vb<=0.0f)&&(vc<=0.0f)) return 0;
2523 if ((va>0.0f)&&(vb>0.0f)&&(vc>0.0f)) return 0;
2524
2525 if ((va<=0.0f)&&(vb>0.0f)&&(vc>0.0f)) return 1; //A
2526 if ((va>0.0f)&&(vb<=0.0f)&&(vc<=0.0f)) return 1;
2527
2528 if ((va>0.0f)&&(vb<=0.0f)&&(vc>0.0f)) return 2;
2529 if ((va<=0.0f)&&(vb>0.0f)&&(vc<=0.0f)) return 2;
2530
2531 if ((va>0.0f)&&(vb>0.0f)&&(vc<=0.0f)) return 3;
2532 if ((va<=0.0f)&&(vb<=0.0f)&&(vc>0.0f)) return 3;
2533
2534 return 4;
2535 }
2536
findContour(QList<V3> & lines,GLfloat value)2537 void FourXle::findContour(QList<V3> &lines, GLfloat value){
2538 GLfloat Va,Vb,Vc,mix1,mix2;
2539 V3 sta,end;
2540 double scope=cScopeBx->value()*cScopeBx->value();
2541 for (int i=0; i<Planpgns.size();i++){
2542 int n=Planpgns.at(i).n;
2543 if (n==3){
2544 Va=Planorte.at(Planpgns.at(i).ii[0]).color-value;
2545 Vb=Planorte.at(Planpgns.at(i).ii[1]).color-value;
2546 Vc=Planorte.at(Planpgns.at(i).ii[2]).color-value;
2547 switch (cutTriangle(Va, Vb, Vc)){
2548 case 1:
2549 mix1=(Va/(Va-Vb));
2550 mix2=(Va/(Va-Vc));
2551 sta=(1.0f-mix1)*Planorte.at(Planpgns.at(i).ii[0]).vertex + (mix1)*Planorte.at(Planpgns.at(i).ii[1]).vertex;
2552 end=(1.0f-mix2)*Planorte.at(Planpgns.at(i).ii[0]).vertex + (mix2)*Planorte.at(Planpgns.at(i).ii[2]).vertex;
2553 if (Distance(sta,end)>1.0) break;
2554 if (scope > 0.0){
2555 if (Distance(sta,aufpunkt)> scope) break;
2556 if (Distance(end,aufpunkt)> scope) break;
2557 }
2558 sta+=pnormal*value*heightScale->value();
2559 end+=pnormal*value*heightScale->value();
2560 lines.append(sta);
2561 lines.append(end);
2562 break;
2563 case 2:
2564 mix1=(Va/(Va-Vb));
2565 mix2=(Vb/(Vb-Vc));
2566 sta=(1.0f-mix1)*Planorte.at(Planpgns.at(i).ii[0]).vertex + (mix1)*Planorte.at(Planpgns.at(i).ii[1]).vertex;
2567 end=(1.0f-mix2)*Planorte.at(Planpgns.at(i).ii[1]).vertex + (mix2)*Planorte.at(Planpgns.at(i).ii[2]).vertex;
2568 if (Distance(sta,end)>1.0) break;
2569 if (cScopeBx->value() > 0.0){
2570 if (Distance(sta,aufpunkt)> scope) break;
2571 if (Distance(end,aufpunkt)> scope) break;
2572 }
2573 sta+=pnormal*value*heightScale->value();
2574 end+=pnormal*value*heightScale->value();
2575 lines.append(sta);
2576 lines.append(end);
2577 break;
2578 case 3:
2579 mix1=(Vc/(Vc-Vb));
2580 mix2=(Va/(Va-Vc));
2581 sta=(1.0f-mix1)*Planorte.at(Planpgns.at(i).ii[2]).vertex + (mix1)*Planorte.at(Planpgns.at(i).ii[1]).vertex;
2582 end=(1.0f-mix2)*Planorte.at(Planpgns.at(i).ii[0]).vertex + (mix2)*Planorte.at(Planpgns.at(i).ii[2]).vertex;
2583 if (Distance(sta,end)>1.0) break;
2584 if (cScopeBx->value() > 0.0){
2585 if (Distance(sta,aufpunkt)> scope) break;
2586 if (Distance(end,aufpunkt)> scope) break;
2587 }
2588 sta+=pnormal*value*heightScale->value();
2589 end+=pnormal*value*heightScale->value();
2590 lines.append(sta);
2591 lines.append(end);
2592 break;
2593 }
2594 }else{
2595 V3 mid=V3(0.0,0.0,0.0);
2596 for (int j=0; j<n;j++) mid+=Planorte.at(Planpgns.at(i).ii[j]).vertex;
2597 mid*=1.0/n;
2598 GLfloat mixColor=0.0f;
2599 for (int j=0; j<n;j++) mixColor+=Planorte.at(Planpgns.at(i).ii[j]).color;
2600 mixColor*=1.0f/n;
2601 Va=mixColor-value;
2602 for (int k=0;k<n;k++){
2603 int kp=(k+1)%n;
2604 Vb=Planorte.at(Planpgns.at(i).ii[k]).color-value;
2605 Vc=Planorte.at(Planpgns.at(i).ii[kp]).color-value;
2606 switch (cutTriangle(Va, Vb, Vc)){
2607 case 1:
2608 mix1=(Va/(Va-Vb));
2609 mix2=(Va/(Va-Vc));
2610 sta=(1.0f-mix1)*mid + (mix1)*Planorte.at(Planpgns.at(i).ii[k]).vertex;
2611 end=(1.0f-mix2)*mid + (mix2)*Planorte.at(Planpgns.at(i).ii[kp]).vertex;
2612 if (Distance(sta,end)>1.0) break;
2613 if (cScopeBx->value() > 0.0){
2614 if (Distance(sta,aufpunkt)> scope) break;
2615 if (Distance(end,aufpunkt)> scope) break;
2616 }
2617 sta+=pnormal*value*heightScale->value();
2618 end+=pnormal*value*heightScale->value();
2619 lines.append(sta);
2620 lines.append(end);
2621 break;
2622 case 2:
2623 mix1=(Va/(Va-Vb));
2624 mix1=(Va/(Va-Vb));
2625 mix2=(Vb/(Vb-Vc));
2626 sta=(1.0f-mix1)*mid + (mix1)*Planorte.at(Planpgns.at(i).ii[k]).vertex;
2627 end=(1.0f-mix2)*Planorte.at(Planpgns.at(i).ii[k]).vertex + (mix2)*Planorte.at(Planpgns.at(i).ii[kp]).vertex;
2628 if (Distance(sta,end)>1.0) break;
2629 if (cScopeBx->value() > 0.0){
2630 if (Distance(sta,aufpunkt)> scope) break;
2631 if (Distance(end,aufpunkt)> scope) break;
2632 }
2633 sta+=pnormal*value*heightScale->value();
2634 end+=pnormal*value*heightScale->value();
2635 lines.append(sta);
2636 lines.append(end);
2637 break;
2638 case 3:
2639 mix1=(Vc/(Vc-Vb));
2640 mix2=(Va/(Va-Vc));
2641 sta=(1.0f-mix1)*Planorte.at(Planpgns.at(i).ii[kp]).vertex + (mix1)*Planorte.at(Planpgns.at(i).ii[k]).vertex;
2642 end=(1.0f-mix2)*mid + (mix2)*Planorte.at(Planpgns.at(i).ii[kp]).vertex;
2643 if (Distance(sta,end)>1.0) break;
2644 if (cScopeBx->value() > 0.0){
2645 if (Distance(sta,aufpunkt)> scope) break;
2646 if (Distance(end,aufpunkt)> scope) break;
2647 }
2648 sta+=pnormal*value*heightScale->value();
2649 end+=pnormal*value*heightScale->value();
2650 lines.append(sta);
2651 lines.append(end);
2652 break;
2653 case 4:
2654 printf("@@~~~~ %g %g %g\n",Va,Vb,Vc);
2655 }
2656 }
2657 }
2658 }
2659 }
2660
CalcPlaneVertex(int ix,int iy,int iz,V3 n,V3 ap,V3 off)2661 void FourXle::CalcPlaneVertex( int ix, int iy, int iz , V3 n, V3 ap, V3 off){
2662 GLfloat vo, vx=0, vy=0, vz=0, mix;
2663 Ort o;
2664 //V3 xpo=V3(1,0,0);
2665 //mole->frac2kart(xpo,xpo);
2666 V3 mdz=0.25*dx+0.25*dy+0.25*dz;
2667 vo = n*(dx*ix + dy*iy + dz*iz-ap); //plane n dot (x - x1)
2668 vx = n*(dx*(ix+1) + dy*iy + dz*iz - ap);
2669 vy = n*(dx*ix + dy*(iy+1) + dz*iz - ap);
2670 vz = n*(dx*ix + dy*iy + dz*(iz+1) - ap);
2671 //V3 fl,mdz=(0.5*dx)+(0.5*dy)+(0.5*dz);
2672 if( ix != n1-1 && Intersect(vo,vx) ){
2673 mix=(vo/(vo-vx));
2674 o.vertex = dx*(mix+ix) + dy*iy + dz*iz +off;//+urs;
2675 /*
2676 mole->kart2frac(o.vertex,o.vertex);
2677 o.vertex+=V3(-0.5,-0.5,-0.5);
2678 fl=V3(floor(o.vertex.x),floor(o.vertex.y),floor(o.vertex.z));
2679 o.vertex+=-1.0*fl;
2680 o.vertex+=V3(0.5,0.5,0.5);
2681 mole->frac2kart(o.vertex,o.vertex);
2682 o.vertex+=-1.0*urs;
2683 // */
2684 o.vertex+=mdz;
2685 //o.vertex+=test3;
2686 //o.vertex+=-1.0*xpo;
2687 o.normal = n;
2688 if (sourceMap->currentIndex()==0)
2689 o.color = datfo_fc[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4]*(1.0-mix)+mix*datfo_fc[((1+ix)%n1)+(iy%n2)*n1+(iz%n3)*n4];
2690 else o.color = datfo[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4]*(1.0-mix)+mix*datfo_fc[((1+ix)%n1)+(iy%n2)*n1+(iz%n3)*n4];
2691 Planorte.append(o);
2692 nodeX[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].index=Planorte.size()-1;
2693 nodeX[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].flag=1;
2694 }else{
2695 nodeX[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].flag = 0;
2696 }
2697 if( iy != n2-1 && Intersect(vo,vy) ){
2698 mix=(vo/(vo-vy));
2699 o.vertex = dx*ix + dy*(mix+iy) + dz*iz+off;//+urs;
2700 /*
2701
2702 mole->kart2frac(o.vertex,o.vertex);
2703 o.vertex+=V3(-0.5,-0.5,-0.5);
2704 fl=V3(floor(o.vertex.x),floor(o.vertex.y),floor(o.vertex.z));
2705 o.vertex+=-1.0*fl;
2706 o.vertex+=V3(0.5,0.5,0.5);
2707 mole->frac2kart(o.vertex,o.vertex);
2708 o.vertex+=-1.0*urs;
2709 // */
2710 o.vertex+=mdz;
2711 //o.vertex+=test3;
2712
2713 //o.vertex+=-1.0*xpo;
2714 o.normal = n;
2715 if (sourceMap->currentIndex()==0)
2716 o.color = datfo_fc[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4]*(1.0-mix)+ mix*datfo_fc[(ix%n1)+((iy+1)%n2)*n1+(iz%n3)*n4];
2717 else o.color = datfo[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4]*(1.0-mix)+ mix*datfo_fc[(ix%n1)+((iy+1)%n2)*n1+(iz%n3)*n4];
2718 Planorte.append(o);
2719 nodeY[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].index=Planorte.size()-1;
2720 nodeY[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].flag=1;
2721 }else{
2722 nodeY[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].flag=0;
2723 }
2724 if( iz != n3-1 && Intersect(vo,vz) ){
2725 mix=(vo/(vo-vz));
2726 o.vertex = dx*ix + dy*iy + (mix+iz)*dz+off;//+urs;
2727 /*
2728 mole->kart2frac(o.vertex,o.vertex);
2729 o.vertex+=V3(-0.5,-0.5,-0.5);
2730 fl=V3(floor(o.vertex.x),floor(o.vertex.y),floor(o.vertex.z));
2731 o.vertex+=-1.0*fl;
2732 o.vertex+=V3(0.5,0.5,0.5);
2733 mole->frac2kart(o.vertex,o.vertex);
2734 o.vertex+=-1.0*urs;
2735 // */
2736
2737 o.vertex+=mdz;
2738 // o.vertex+=test3;
2739 //o.vertex+=-1.0*xpo;
2740 o.normal = n;
2741 if (sourceMap->currentIndex()==0)
2742 o.color = datfo_fc[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4]*(1.0-mix)+mix*datfo_fc[(ix%n1)+(iy%n2)*n1+((iz+1)%n3)*n4];
2743 else o.color = datfo[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4]*(1.0-mix)+mix*datfo_fc[(ix%n1)+(iy%n2)*n1+((iz+1)%n3)*n4];
2744 Planorte.append(o);
2745 nodeZ[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].index = Planorte.size()-1;
2746 nodeZ[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].flag=1;
2747 }else{
2748 nodeZ[(ix%n1)+(iy%n2)*n1+(iz%n3)*n4].flag = 0;
2749 }
2750 }
2751
makeFacesp(int nn,Node poly[])2752 void FourXle::makeFacesp(int nn, Node poly[] ){
2753 // /*
2754 int n=0;
2755 int ly[13];
2756 for (int j=0;j<nn;j++)//zu nahe beieinander liegende verts ignorieren
2757 if (Distance(Planorte.at(poly[j].index).vertex,Planorte.at(poly[(j+1)%nn].index).vertex)>0.00)
2758 ly[n++]=poly[j].index;
2759 else{
2760 Planorte[poly[j].index].vertex=
2761 Planorte[poly[(j+1)%nn].index].vertex=
2762 (Planorte.at(poly[j].index).vertex+Planorte.at(poly[(j+1)%nn].index).vertex)*0.5f;
2763 Planorte[poly[j].index].color=
2764 Planorte[poly[(j+1)%nn].index].color=
2765 (Planorte.at(poly[j].index).color+Planorte.at(poly[(j+1)%nn].index).color)*0.5f;
2766 ly[n++]=poly[(j+1)%nn].index;
2767 j++;
2768 }
2769 GLfloat mid_col = 0.0;
2770 V3 mid_ver = V3(0,0,0);
2771 V3 mid_nor = V3(0,0,0);
2772 for(int i=0; i<n; i++ ){
2773 mid_ver += Planorte.at(ly[i]).vertex;
2774 mid_nor += Planorte.at(ly[i]).normal;
2775 mid_col += Planorte.at(ly[i]).color;
2776 }
2777 mid_ver *= (1.0/n);
2778 Normalize(mid_nor);
2779 mid_col *= (1.0/n);
2780 if (n<3) return; //weniger als 3 verts -> nichts zu tun
2781 Polygn neupoly;
2782 neupoly.vertex=mid_ver;
2783 neupoly.n=n;
2784 //if (Planorte.at(poly[0].index).direct*iso_level<0.0){//drehsinn der polygone abfragen* /
2785 for(int i=0; i<n; i++ ){
2786 neupoly.ii[i]=ly[i];
2787 }
2788 //}
2789 /*else{
2790 for(int i=n-1,k=0; i>=0; i-- ){
2791 neupoly.ii[i]=ly[k++];
2792 }
2793 }//else*/
2794 Planpgns.append(neupoly);
2795 // */
2796 }
2797
MakeElementp(int ix,int iy,int iz,int s1,int s2)2798 void FourXle::MakeElementp( int ix, int iy, int iz ,int s1, int s2) {//das ist der Teil des japanischen Programms den ich nicht verstehe.
2799 //Hauptsache fuktioniert.
2800 // /*
2801 static int conn[12][2][4] = { //char->int wegen warning g++ >3.0
2802 {{ 0, 1, 7, 6}, { 0, 2, 8, 3}}, // 0
2803 {{ 1, 2, 5, 4}, { 1, 0, 6, 7}}, // 1
2804 {{ 2, 0, 3, 8}, { 2, 1, 4, 5}}, // 2
2805 {{ 3, 8, 2, 0}, { 3, 4,10, 9}}, // 3
2806 {{ 4, 3, 9,10}, { 4, 5, 2, 1}}, // 4
2807 {{ 5, 4, 1, 2}, { 5, 6, 9,11}}, // 5
2808 {{ 6, 5,11, 9}, { 6, 7, 1, 0}}, // 6
2809 {{ 7, 6, 0, 1}, { 7, 8,11,10}}, // 7
2810 {{ 8, 7,10,11}, { 8, 3, 0, 2}}, // 8
2811 {{ 9,10, 4, 3}, { 9,11, 5, 6}}, // 9
2812 {{10,11, 8, 7}, {10, 9, 3, 4}}, // 10
2813 {{11, 9, 6, 5}, {11,10, 7, 8}} // 11
2814 };
2815 static Node node[12];
2816 static Node polygon[12];
2817
2818 node[ 0] = nodeX[(ix+iy*s1+iz*s2) %n5];// 000x
2819 node[ 1] = nodeY[(ix+iy*s1+iz*s2) %n5];// 000y
2820 node[ 2] = nodeZ[(ix+iy*s1+iz*s2) %n5];// 000z
2821 node[ 3] = nodeX[(ix+iy*s1+((iz+1)%n3)*s2) %n5];// 001y
2822 node[ 4] = nodeY[(ix+iy*s1+((iz+1)%n3)*s2) %n5];// 001z
2823 node[ 5] = nodeZ[(ix+((iy+1)%n2)*s1+iz*s2) %n5];// 010x
2824 node[ 6] = nodeX[(ix+((iy+1)%n2)*s1+iz*s2) %n5];// 010y
2825 node[ 7] = nodeY[(((1+ix)%n1)+iy*s1+iz*s2) %n5];// 100y
2826 node[ 8] = nodeZ[(((1+ix)%n1)+iy*s1+iz*s2) %n5];// 100z
2827 node[ 9] = nodeX[(ix+((iy+1)%n2)*s1+((iz+1)%n3)*s2)%n5];// 011x
2828 node[10] = nodeY[(((ix+1)%n1)+iy*s1+((iz+1)%n3)*s2)%n5];// 101y
2829 node[11] = nodeZ[(((ix+1)%n1)+((iy+1)%n2)*s1+iz*s2)%n5];// 110z
2830
2831 if (((char)node[0]+node[1]+node[2]+node[3]+node[4]+node[5]
2832 +node[6]+node[7]+node[8]+node[9]+node[10]+node[11])==0) return;
2833 for( int is=0; is<12; is++ ) {
2834 if( !node[is] ) continue;
2835
2836 int n=0, i=is, m=0;//,ai=i;
2837 do {
2838 polygon[n++]= node[i];
2839 int sol = IndexSelectedP(node[conn[i][m][0]],node[conn[i][m][1]],
2840 node[conn[i][m][2]],node[conn[i][m][3]]);
2841 // ai=i;
2842 i = conn[i][m][sol];
2843
2844 if( sol == 2 ) m ^= 1;
2845 node[i].flag = 0;
2846 }while( (i!=is)&&(n<11) );
2847 {
2848 makeFacesp( n, polygon );
2849 }
2850 }
2851 // */
2852 }
2853
makePlane(QList<V3> & lines,int a1,int a2,int a3)2854 void FourXle::makePlane(QList<V3> &lines,int a1, int a2, int a3) {
2855 /*
2856 farbe[0][0]=0.6;
2857 farbe[0][1]=0.0;
2858 farbe[0][2]=0.0;
2859 farbe[0][3]=1.0;
2860
2861 farbe[1][0]=0.6;
2862 farbe[1][1]=0.0;
2863 farbe[1][2]=0.0;
2864 farbe[1][3]=1.0;
2865
2866 farbe[2][0]=0.6;
2867 farbe[2][1]=0.0;
2868 farbe[2][2]=0.0;
2869 farbe[2][3]=1.0;
2870
2871 farbe[3][0]=0.0;
2872 farbe[3][1]=0.0;
2873 farbe[3][2]=0.0;
2874 farbe[3][3]=1.0;
2875
2876 farbe[4][0]=0.0;
2877 farbe[4][1]=0.0;
2878 farbe[4][2]=1.0;
2879 farbe[4][3]=1.0;
2880
2881 farbe[5][0]=0.0;
2882 farbe[5][1]=0.0;
2883 farbe[5][2]=1.0;
2884 farbe[5][3]=1.0;
2885
2886 farbe[6][0]=0.0;
2887 farbe[6][1]=0.0;
2888 farbe[6][2]=1.0;
2889 farbe[6][3]=1.0;
2890
2891 Farben=7;
2892 */
2893
2894 // fprintf(stderr, "makePlane %d %d %d %d %d\n",n1,n2,n3,n4,n5);
2895 //orig=V3(0,0,0);
2896 // n4=n1*n2;
2897 // data.clear();
2898
2899 if (( nodeX =(Node*)malloc(sizeof(Node)*n5))==nullptr) {
2900 qDebug()<<"Less Memory(X) ";
2901 exit(1);
2902 }
2903 if (( nodeY =(Node*)malloc(sizeof(Node)*n5))==nullptr) {
2904 qDebug()<<"Less Memory(Y)!!";
2905 exit(1);
2906 }
2907 if (( nodeZ =(Node*)malloc(sizeof(Node)*n5))==nullptr) {
2908 qDebug()<<"Less Memory(Z)!!";
2909 exit(1);
2910 }
2911 // printf("UIUI %p %p %p\n",nodeX,nodeY,nodeZ);
2912 for (int ix=0; ix<n5; ix++){
2913 nodeX[ix].index=0;
2914 nodeX[ix].flag=0;
2915 nodeY[ix].index=0;
2916 nodeY[ix].flag=0;
2917 nodeZ[ix].index=0;
2918 nodez[ix].flag=0;
2919 }
2920 // printf("__\n");
2921 Planorte.clear();
2922 Planpgns.clear();
2923 lines.clear();
2924 chgl->contval.clear();
2925
2926 //test3= ((n1-1)/-2.0) * dx + ((n2-1)/-2.0) * dy + ((n3-1)/-2.0) * dz + orig;
2927 //test3=orig;
2928 // printf("Atomindices: %d %d %d < %d\n",a1,a2,a3,mole->showatoms.size());
2929 if ((a1>=mole->showatoms.size())||(a2>=mole->showatoms.size())||(a3>=mole->showatoms.size())) return;
2930 if ((a1<0)||(a2<0)||(a3<0)) return;
2931 if (mole->showatoms[a1].pos==mole->showatoms[a2].pos){
2932 mole->frac2kart(mole->showatoms[a1].frac,mole->showatoms[a1].pos);
2933 mole->frac2kart(mole->showatoms[a2].frac,mole->showatoms[a2].pos);
2934 mole->frac2kart(mole->showatoms[a3].frac,mole->showatoms[a3].pos);
2935 }
2936 V3 a1v=V3(mole->showatoms.at(a1).pos.x,mole->showatoms.at(a1).pos.y,mole->showatoms.at(a1).pos.z);
2937 V3 a2v=V3(mole->showatoms.at(a2).pos.x,mole->showatoms.at(a2).pos.y,mole->showatoms.at(a2).pos.z);
2938 V3 a3v=V3(mole->showatoms.at(a3).pos.x,mole->showatoms.at(a3).pos.y,mole->showatoms.at(a3).pos.z);
2939 pnormal = Normalize((a2v - a1v) % (a3v - a1v));
2940 switch (centerIsOn->currentIndex()){
2941 case 0: aufpunkt = V3(mole->showatoms.at(a1).pos.x,mole->showatoms.at(a1).pos.y,mole->showatoms.at(a1).pos.z);break;
2942 case 1: aufpunkt = V3(mole->showatoms.at(a2).pos.x,mole->showatoms.at(a2).pos.y,mole->showatoms.at(a2).pos.z);break;
2943 case 2: aufpunkt = V3(mole->showatoms.at(a3).pos.x,mole->showatoms.at(a3).pos.y,mole->showatoms.at(a3).pos.z);break;
2944 case 3: {
2945 aufpunkt = (1.0/3.0)*
2946 (V3(mole->showatoms.at(a1).pos.x,mole->showatoms.at(a1).pos.y,mole->showatoms.at(a1).pos.z)+
2947 V3(mole->showatoms.at(a2).pos.x,mole->showatoms.at(a2).pos.y,mole->showatoms.at(a2).pos.z)+
2948 V3(mole->showatoms.at(a3).pos.x,mole->showatoms.at(a3).pos.y,mole->showatoms.at(a3).pos.z));
2949 break;
2950 }
2951
2952 }
2953 // aufpunkt = V3(mole->showatoms.at(a1).kart.x,mole->showatoms.at(a1).kart.y,mole->showatoms.at(a1).kart.z);
2954 // printf("%g %g %g %g %g %g %g\n",Norm(pnormal),pnormal.x,pnormal.y,pnormal.z,aufpunkt.x,aufpunkt.y,aufpunkt.z);
2955 double angle = winkel(pnormal, V3(0,0,1));
2956 V3 ax= Normalize(pnormal % V3(0,0,1));
2957 double cp = cos(angle), cp1 = 1.0 - cos(angle), sp = sin(angle);
2958 double X=ax.x,Y=ax.y,Z=ax.z;
2959 Matrix ROT = Matrix(
2960 cp+cp1*X*X, -Z*sp+cp1*X*Y, Y*sp+cp1*X*Z,
2961 Z*sp+cp1*X*Y, cp+cp1*Y*Y, -X*sp+cp1*Y*Z,
2962 -Y*sp+cp1*X*Z, X*sp+cp1*Y*Z, cp+cp1*Z*Z);
2963 V3 bv=ROT*Normalize(V3(a2v.x-a1v.x,a2v.y-a1v.y,a2v.z-a1v.z));
2964 //printf("BV %g %g %g\n",bv.x,bv.y,bv.z);
2965 double bww=winkel(bv,V3( 1,0,0));
2966 cp=cos(bww);
2967 sp=sin(bww);
2968 //printf("%g %g %g\n",bww*180.0/M_PI,cp,sp);
2969 Matrix R2 = Matrix(
2970 cp, sp, 0,
2971 -sp, cp, 0,
2972 0, 0, 1);//*/
2973 bv=R2*bv;
2974 bww=winkel(bv,V3( 1,0,0));
2975 if (fabs(bww)>1) {
2976 R2 = Matrix(
2977 cp,-sp, 0,
2978 sp, cp, 0,
2979 0, 0, 1);//*/
2980
2981 }
2982 //printf("BV %g %g %g %g\n",bv.x,bv.y,bv.z,bww);
2983 V3 vt = ROT * V3(pnormal.x, pnormal.y, pnormal.z);
2984 vt = R2*vt;
2985 //printf("ang %g %g %g \n%12.6f %12.6f %12.6f \n%g %g\n%g\n",angle/M_PI*180.0, pnormal*ax,ax*V3(0,0,1), vt.x, vt.y, vt.z,determinant(ROT),Norm(ax),winkel(a1v,V3(0,1,0)));
2986 // fprintf(stderr, "makePlane\n");
2987 int qii=0;V3 ap;
2988 int aa=-1,ee=2;
2989 if (!extendUnit->isChecked()){aa++;ee--;}
2990 for (int iix=aa; iix<ee; iix++){
2991 for (int iiy=aa; iiy<ee; iiy++){
2992 for (int iiz=aa; iiz<ee; iiz++){
2993 V3 propo=V3(iix,iiy,iiz)+V3(0.5,0.5,0.5);
2994 V3 pp=V3(iix,iiy,iiz);
2995 mole->frac2kart(propo,propo);
2996 mole->frac2kart(pp,pp);
2997 double abst=pnormal*(propo-aufpunkt);
2998 V3 da=propo-abst*pnormal;
2999 mole->kart2frac(da,da);
3000 qii++;
3001 // printf("PLANE %2d %2d %2d %12g %12g %12g %g %12g %12g %12g\n",iix,iiy,iiz,da.x,da.y,da.z,abst,da.x-iix,da.y-iiy,da.z-iiz);
3002 ap=da-V3(iix,iiy,iiz);
3003 mole->frac2kart(ap,ap);
3004 //printf("H%d' 1 %8.5f %8.5f %8.5f 11.0 0.05 %5.3f \n",qii, da.x, da.y, da.z,abst);
3005 // if ((da.x-iix<1.0)&&(da.x-iix>0.0)&&(da.y-iiy<1.0)&&(da.y-iiy>0.0)&&(da.z-iiz<1.0)&&(da.z-iiz>0.0)){
3006 for( int ix=0; ix<=n1; ix++ ){
3007 for( int iy=0; iy<=n2; iy++ ){
3008 for( int iz=0; iz<=n3; iz++ ){
3009 CalcPlaneVertex(ix,iy,iz,pnormal,ap,pp);//
3010 //CalcPlaneVertex(ix,iy,iz,pnormal,aufpunkt,V3(0,0,0));
3011 }
3012 }
3013 }
3014 //printf("%3d %3d %3d Planorte.size = %d\n",iix,iiy,iiz,Planorte.size());
3015 //fprintf(stderr, "makePlane\n");
3016 for( int ix=0; ix<n1-1; ix++ )
3017 for( int iy=0; iy<n2-1; iy++ )
3018 for( int iz=0; iz<n3-1; iz++ )
3019 MakeElementp(ix,iy,iz,n1,n4);
3020 //}
3021 }
3022 }
3023 }
3024 double minP=1e37,maxP=-1e38;
3025 for (int i=0; i<Planorte.size(); i++) {
3026 minP=fmin(Planorte.at(i).color,minP);
3027 maxP=fmax(Planorte.at(i).color,maxP);
3028 }
3029 //fprintf(stderr, "makePlane Planpgns.size = %d Planorte.size = %d minP = %g maxP = %g\n",Planpgns.size(),Planorte.size(),minP,maxP);
3030 min= 1e37;
3031 max=-1e37;
3032 QStringList numbers = QString(contourValueEdit->text()).split(" ",skipEmptyParts);
3033 int z=0;
3034 // printf("ok\n");
3035 if (!numbers.isEmpty()) chgl->contval[0]=numbers.at(0).toFloat();
3036 for (int i = 0; i < numbers.size(); i++){
3037 float contour = numbers.at(i).toFloat();
3038 findContour(lines,contour);
3039 chgl->contval[lines.size()]=contour;
3040 if ((lines.size()-z)>0){
3041 min=fmin(contour,min);
3042 max=fmax(contour,max);
3043 }
3044 z=lines.size();
3045 }
3046 chgl->contMax=max;
3047 chgl->contMin=min;
3048 // printf("ok %f %f\n",min,max);
3049 //legende();
3050 free(nodeX);
3051 nodeX=nullptr;
3052 free(nodeY);
3053 nodeY=nullptr;
3054 free(nodeZ);
3055 nodeZ=nullptr;
3056 Planpgns.clear();
3057 Planorte.clear();
3058 //printf("ok\n");
3059 QList<V3> lin2d;
3060 int aspect1=3,aspect2=4;//todo make aspect ratio changable
3061 switch (aspectRatios->currentIndex()) {
3062 case 0: aspect1=3,aspect2=4; break;
3063 case 1: aspect1=9,aspect2=16; break;
3064 case 2: aspect1=1,aspect2=1; break;
3065 }
3066 int width=1024*aspect1/aspect2;
3067 double miX = 9.e30;
3068 double maX = -9.e30;
3069 double miY = 9.e30;
3070 double maY = -9.e30;
3071 V3 v1;//,v2=V3(99,99,99);
3072 for (int i=0; i<lines.size(); i++){
3073 v1=R2*(ROT * V3(lines.at(i).x-aufpunkt.x,lines.at(i).y-aufpunkt.y,lines.at(i).z-aufpunkt.z));
3074 //if (Distance(v1,v2)>0.0)
3075 lin2d.append(v1);
3076 miX=fmin(v1.x,miX);
3077 maX=fmax(v1.x,maX);
3078 miY=fmin(v1.y,miY);
3079 maY=fmax(v1.y,maY);
3080
3081
3082 //v2=v1;
3083 }
3084 //printf("!%d == %d! %f < X < %f, %f < Y < %f \n",lin2d.size(),lines.size(),miX,maX,miY,maY);
3085 if (cScopeBx->value() > 0.0){
3086 double XX=cScopeBx->value();
3087 double YY=(XX*aspect1)/aspect2;//todo make aspect ratio changable
3088 miX=-XX;
3089 miY=-YY;
3090 maX= XX;
3091 maY= YY;
3092 }
3093 //printf("!%d == %d! %f < X < %f, %f < Y < %f \n",lin2d.size(),lines.size(),miX,maX,miY,maY);
3094 if (contEPSFile->text().isEmpty()){
3095 // qDebug()<<"contEPSFile->text() is Empty()";
3096 QString contourDescription1;
3097 for (int i=0; i<lin2d.size()/2; i++){
3098 //(wrt-min)/(max-min)
3099 if (chgl->contval.contains(i*2)) {
3100 if (chgl->contval.value(i*2)>0.00001f) {
3101 contourDescription1.append(QString("<font color=blue>blue contours at %1 eÅ<sup>-3</sup></font> <br>\n").arg(chgl->contval.value(i*2)));
3102 }
3103 else if (chgl->contval.value(i*2)<-0.00001f){
3104 contourDescription1.append(QString("<font color=red>red contours at %1 eÅ<sup>-3</sup></font> <br>\n").arg(chgl->contval.value(i*2)));
3105 }
3106 else {
3107 contourDescription1.append(QString("<font color=black>black contours at 0 eÅ<sup>-3</sup></font> <br>\n"));
3108 }
3109 }
3110 }
3111 emit bigmessage(QString("%1<br>").arg(contourDescription1));
3112 return;
3113 }
3114
3115 // ||(!QFile::exists (contEPSFile->text())))
3116 //qDebug()<<QDate::currentDate().toString(Qt::ISODate);
3117 FILE *f=fopen(contEPSFile->text().toLocal8Bit(),"wt");
3118 fprintf(f,"%s%d\n","%!PS-Adobe-3.0 EPSF-3.0\n%%BoundingBox: 0 0 1024 ",width);
3119 fprintf(f,"%s",
3120 "%%Title: %s Contour Plot\n"
3121 "%%Creator: %s by Christian B. Hubschle\n",PROGRAM_NAME,PROGRAM_NAME);
3122 fprintf(f,"%%%%CreationDate: %s\n",QDate::currentDate().toString(Qt::ISODate).toStdString().c_str());
3123
3124 fprintf(f,"%s",
3125 //"%%Pages: 1\n"
3126 "%%EndComments\n"
3127 "%%BeginProlog\n"
3128 //"/l { newpath moveto lineto stroke } bind def\n"
3129 //"/nm { newpath moveto } bind def\n"
3130 "/cp {closepath} bind def\n"
3131 "/ef {eofill} bind def\n"
3132 "/gr {grestore} bind def\n"
3133 "/gs {gsave} bind def\n"
3134 "/sa {save} bind def\n"
3135 "/rs {restore} bind def\n"
3136 "/l {lineto} bind def\n"
3137 "/m {moveto} bind def\n"
3138 "/rm {rmoveto} bind def\n"
3139 "/n {newpath} bind def\n"
3140 "/s {stroke} bind def\n"
3141 "/sh {show} bind def\n"
3142 "/slc {setlinecap} bind def\n"
3143 "/slj {setlinejoin} bind def\n"
3144 "/slw {setlinewidth} bind def\n"
3145 "/srgb {setrgbcolor} bind def\n"
3146 "/rot {rotate} bind def\n"
3147 "/sc {scale} bind def\n"
3148 "/sd {setdash} bind def\n"
3149 "/ff {findfont} bind def\n"
3150 "/sf {setfont} bind def\n"
3151 "/scf {scalefont} bind def\n"
3152 "/sw {stringwidth} bind def\n"
3153 //"%%EndSetup\n"
3154 "%%EndProlog\n"
3155 //"%%Page: 1 1\n"
3156 //"gsave\n"
3157 //"20 20 scale\n"
3158 //"10 10 translate\n"
3159 "save\n0 slc 0 slj\n/Helvetica ff 18 scf sf\n");
3160 fprintf(f,"n 0 0 m 1024 0 l 1024 %d l 0 %d l cp clip s\n",width,width);
3161 fprintf(f,"n 2 2 m 1023 2 l 1023 %d l 2 %d l cp s \n",width-2,width-2);
3162 //"0 setgray\n"
3163
3164
3165 //Atoms
3166 for (int i=0; i<mole->showatoms.size(); i++){
3167 V3 ato=R2*(ROT*(mole->showatoms.at(i).pos-V3(aufpunkt.x,aufpunkt.y,aufpunkt.z)));
3168 if ((ato.x>maX)||(ato.x<miX)||(ato.y>maY)||(ato.y<miY)) continue;
3169 //Acol[xdinp[j].OrdZahl]
3170 if (fabs(ato.z)<0.2){
3171 double rad=(mole->showatoms.at(i).an>-1)?mole->arad[mole->showatoms.at(i).an]:0.15;
3172 QColor col=mole->AtomColor[mole->showatoms.at(i).an];
3173 fprintf(f,"n %G %G %G 0 360 arc cp gs %g %g %g srgb fill s gr%% %g %s\n",
3174 (ato.x-miX)/(maX-miX)*1024, (ato.y-miY)/(maY-miY)*width, rad*50,
3175 col.redF(),
3176 col.greenF(),
3177 col.blueF(),ato.z,
3178 mole->showatoms.at(i).Label.toStdString().c_str());
3179 }
3180 // 4 5 3 0 360 arc closepath
3181 }
3182 //Bonds
3183 for (int k=0;k<mole->showbonds.size();k++){
3184 V3 anf=R2*(ROT*(mole->showbonds[k].ato1->pos-V3(aufpunkt.x,aufpunkt.y,aufpunkt.z)));
3185 V3 end=R2*(ROT*(mole->showbonds[k].ato2->pos-V3(aufpunkt.x,aufpunkt.y,aufpunkt.z)));
3186 if ((anf.x>maX)||(anf.x<miX)||(anf.y>maY)||(anf.y<miY)) continue;
3187 if ((end.x>maX)||(end.x<miX)||(end.y>maY)||(end.y<miY)) continue;
3188 if ((fabs(anf.z)<0.2)&&(fabs(end.z)<0.2))
3189 fprintf(f,"n %G %G m %G %G l cp s\n",
3190 (anf.x-miX)/(maX-miX)*1024, (anf.y-miY)/(maY-miY)*width,
3191 (end.x-miX)/(maX-miX)*1024, (end.y-miY)/(maY-miY)*width);
3192
3193 }
3194
3195 // V3 v1,v2;
3196 //qDebug()<<contval;
3197 float red=0.0f,green=0.0f,blue=0.0f;
3198 QString contourDescription;
3199 for (int i=0; i<lin2d.size()/2; i++){
3200 //(wrt-min)/(max-min)
3201 if (chgl->contval.contains(i*2)) {
3202 if (chgl->contval.value(i*2)>0.00001f) {
3203 red=0.0f;green=0.0f;blue=1.0f;
3204 contourDescription.append(QString("% blue@%1, %2 lines\n").arg(chgl->contval.value(i*2)).arg(i));
3205 // printf("blau %f %d\n",contval.value(i*2),i);
3206 }
3207 else if (chgl->contval.value(i*2)<-0.00001f){
3208 red=1.0f;green=0.0f;blue=0.0f;
3209 contourDescription.append(QString("% red@%1, %2 lines\n").arg(chgl->contval.value(i*2)).arg(i));
3210 // printf("rot %f %d\n",contval.value(i),i);
3211 }
3212 else {
3213 red=0.0f;green=0.0f;blue=0.0f;
3214 contourDescription.append(QString("% black@%1, %2 lines\n").arg(chgl->contval.value(i*2)).arg(i));
3215 // printf("schwarz %f %d\n",chgl->contval.value(i*2),i);
3216 }
3217 }
3218 if ((lin2d.at(2*i ).x>maX)||(lin2d.at(2*i ).x<miX)||(lin2d.at(2*i ).y>maY)||(lin2d.at(2*i ).y<miY)) continue;
3219 if ((lin2d.at(2*i+1).x>maX)||(lin2d.at(2*i+1).x<miX)||(lin2d.at(2*i+1).y>maY)||(lin2d.at(2*i+1).y<miY)) continue;
3220 fprintf(f,"n %G %G m %G %G l cp gs %g %g %g srgb 0.5 slw s gr\n",
3221 (lin2d.at(2*i ).x-miX)/(maX-miX)*1024, (lin2d.at(i*2 ).y-miY)/(maY-miY)*width,
3222 (lin2d.at(2*i+1).x-miX)/(maX-miX)*1024, (lin2d.at(i*2+1).y-miY)/(maY-miY)*width,red,green,blue);
3223 }
3224 for (int i=0; i<mole->showatoms.size(); i++){
3225 V3 ato=R2*(ROT*(mole->showatoms.at(i).pos-V3(aufpunkt.x,aufpunkt.y,aufpunkt.z)));
3226 if ((ato.x>maX)||(ato.x<miX)||(ato.y>maY)||(ato.y<miY)) continue;
3227 if (fabs(ato.z)<0.2){
3228 double rad=(mole->showatoms.at(i).an>-1)?mole->arad[mole->showatoms.at(i).an]:0.15;
3229 fprintf(f,"gs n %G %G m (%s) true charpath 1 slw 1 setgray s gr\n",(ato.x-miX)/(maX-miX)*1024+rad*50, (ato.y-miY)/(maY-miY)*width,mole->showatoms.at(i).Label.toStdString().c_str());
3230 fprintf(f,"gs %G %G m (%s) sh gr\n",(ato.x-miX)/(maX-miX)*1024+rad*50, (ato.y-miY)/(maY-miY)*width,mole->showatoms.at(i).Label.toStdString().c_str());
3231
3232 }
3233 // 4 5 3 0 360 arc closepath
3234 }
3235 fprintf(f,"%s",
3236 "restore\n"
3237 "showpage\n"
3238 "%%Trailer\n"
3239 "%%EOF\n"
3240 );
3241 fprintf(f,"%% plane: %s %s %s\n%% contours: %s\n",
3242 mole->showatoms.at(a1).Label.toStdString().c_str(),
3243 mole->showatoms.at(a2).Label.toStdString().c_str(),
3244 mole->showatoms.at(a3).Label.toStdString().c_str(),
3245 contourDescription.toStdString().c_str());
3246 //printf("->%d\n",fclose(f));
3247 //printf("ok3074\n");
3248 QString contourDescription1;
3249 for (int i=0; i<lin2d.size()/2; i++){
3250 //(wrt-min)/(max-min)
3251 if (chgl->contval.contains(i*2)) {
3252 if (chgl->contval.value(i*2)>0.00001f) {
3253 contourDescription1.append(QString("<font color=blue>blue contours at %1 eÅ<sup>-3</sup></font> <br>\n").arg(chgl->contval.value(i*2)));
3254 }
3255 else if (chgl->contval.value(i*2)<-0.00001f){
3256 contourDescription1.append(QString("<font color=red>red contours at %1 eÅ<sup>-3</sup></font> <br>\n").arg(chgl->contval.value(i*2)));
3257 }
3258 else {
3259 contourDescription1.append(QString("<font color=black>black contours at 0 eÅ<sup>-3</sup></font> <br>\n"));
3260 }
3261 }
3262 }
3263 emit bigmessage(QString("%1<br>").arg(contourDescription1));
3264
3265 }
3266
3267
3268
3269
